Major Histocompatibility Complex Class II (DRB3) Genetic Diversity in Spanish Morucha and Colombian Normande Cattle Compared to Taurine and Zebu Populations

Bovine leukocyte antigens (BoLA) have been used as disease markers and immunological traits in cattle due to their primary role in pathogen recognition by the immune system. A higher MHC allele diversity in a population will allow presenting a broader peptide repertoire. However, loss of overall diversity due to domestication process can decrease a population’s peptide repertoire. Within the context of zebu and taurine cattle populations, BoLA-DRB3 genetic diversity in Spanish Morucha and Colombian Normande cattle was analyzed and an approach to estimate functional diversity was performed. Sequence-based typing was used for identifying 29, 23, 27, and 28 alleles in Spanish Morucha, Nariño-, Boyacá-, and Cundinamarca-Normande cattle, respectively. These breeds had remarkably low heterozygosity levels and the Hardy-Weinberg principle revealed significant heterozygote deficiency. FST and DA genetic distance showed that Colombian Normande populations had greater variability than other phenotypically homogeneous breeds, such as Holstein. It was also found that Spanish Morucha cattle were strongly differentiated from other cattle breeds. Spanish Morucha had greater divergence in the peptide-binding region regarding other cattle breeds. However, peptide-binding region covariation indicated that the potential peptide repertoire seemed equivalent among cattle breeds. Despite the genetic divergence observed, the extent of the potential peptide repertoire in the cattle populations studied appears to be similar and thus their pathogen recognition potential should be equivalent, suggesting that functional diversity might persist in the face of bottlenecks imposed by domestication and breeding.


INTRODUCTION
The major histocompatibility complex (MHC) is a primary component of the adaptive immune system that offers a unique alternative for addressing both immunological and evolutionary biological issues. The MHC consists of a group of loci in jawed vertebrates encoding molecules which are fundamental for regulating the immune response (Klein et al., 2007). Class I and class II genes comprise one subset of these loci encoding the cell surface glycoproteins necessary for Tlymphocytes to recognize antigenic peptides on cell surfaces (Rock et al., 2016). While class I molecules are expressed by all nucleated cells and present peptides from intracellular proteins to CD8 + T-cells, class II molecules are expressed by professional antigen-presenting cells and present peptides derived from extracellular proteins to CD4 + T-cells, thereby playing a key role in defence against pathogens (Klein et al., 2007;Neefjes et al., 2011).
The MHC in cattle (also known as bovine leukocyte antigen-BoLA) is located on chromosome 23 and has the overall structure characterizing other mammals' MHC (Takeshima and Aida, 2006). DR and DQ genes encode the molecules within the BoLA class II region that bind the peptides which will be presented to T-lymphocytes, representing the main class II restriction elements for CD4 + T-cells. The BoLA-DQ region comprises DQA and DQB loci, which may vary in number depending on the haplotype; this enables additional diversity by means of intrahaplotype and interhaplotype pairing of DQA and DQB molecules (Andersson and Rask, 1988;Glass et al., 2000;Norimine and Brown, 2005). BoLA-DR consists of the monomorphic BoLA-DRA locus and three DRB loci, of which the BoLA-DRB3 gene is the only one known to be fully functional (Burke et al., 1991).
BoLA-DRB3 is the most polymorphic bovine MHC gene; 136 different alleles have been reported to date (Maccari et al., 2017). Such polymorphisms are mainly located in the b1 domain's peptide-binding region (PBR), which is encoded by exon 2 and has been used for defining BoLA-DRB3 alleles (Sigurdardottir et al., 1991). Because the amino acids (aa) forming the PBR determine MHC-presented peptides' binding affinity, different alleles will bind a different repertoire of peptides, thus influencing immune response variability (Baxter et al., 2009;van Deutekom and Kesmir, 2015). Similarities between BoLA-DRB3 PBR are used as criteria for developing in silico panspecific methods for predicting MHC-peptide-binding affinity (Nielsen et al., 2008;Hoof et al., 2009). BoLA-DRB3 divergence can be used for estimating the size of peptide-binding repertoire. Thus, individuals or populations having highly divergent BoLA-DRB3 alleles will have a broader peptide-binding repertoire than those having very similar alleles (Klein et al., 2007;Lenz et al., 2013). Furthermore, peptide-binding differences among populations could be estimated by correlating variation patterns among their BoLA-DRB3 allele repertoires.
MHC diversity is of great interest for breeders, population geneticists and evolutionary biologists, considering that different degrees and MHC variability patterns reflect evolutionary processes such as adaptation, selection (natural, sexual or artificial) and drift within and between populations (Goszczynski et al., 2014;Takeshima et al., 2014). Several studies have shown that decreased MHC variability might be caused by population bottlenecks (Bollmer et al., 2007;Bollmer et al., 2011;Mason et al., 2011;Zhang et al., 2016). On the contrary, a high level of diversity could be maintained by balancing selection driven by pathogens or other mechanisms due to MHC function regarding pathogen recognition, despite extreme population bottlenecks (Edwards and Hedrick, 1998;Garrigan and Hedrick, 2001;Aguilar et al., 2004;Borg et al., 2011;Moutou et al., 2013;Newhouse and Balakrishnan, 2015). Although functional polymorphism may be narrower than genetic polymorphism in pigs (Moutou et al., 2013), most studies on domestic animals' MHC have focused on genetic polymorphism; however, it remains unknown whether bottlenecks associated with domestication and breeding could have reduced MHC functional diversity.
Despite advances in characterizing BoLA-DRB3 genetic diversity in cattle, information about allele frequency is not available for some breeds (only for a few out of over 800 breeds of cattle recognized worldwide via sequence-based typing) (Takeshima et al., 2018). The few MHC alleles in domestic animals contrasts with those of other species (Nino-Vasquez et al., 2000;Suarez et al., 2006;Lopez et al., 2014;Maccari et al., 2017;Shiina et al., 2017). However, a smaller amount of alleles does not necessarily mean reduced functionality, in fact two species can have very similar functional repertoire size despite varying regarding the amount of alleles (Suarez et al., 2017).
This study has thus been aimed at determining whether domestic cattle breeding has decreased actual MHC variability. MHC diversity of two new breeds has been characterized. The Normande breed was chosen out of convenience. Spanish Morucha, a breed having relatively low human intervention, has been used here as an approach to what could be considered a natural cattle population. BoLA-DRB3 genetic diversity, structure, differentiation, and selective pressure in Colombian Normande and Spanish Morucha breeds was characterized and compared to worldwide taurine and zebu breeds. Each cattle breed's BoLA-DRB3 PBR variability and covariation was then analyzed as an approach to compare potential peptide-binding repertoire size and indication of functional diversity. Genetic divergence patterns based on population genetics and PBR variability analysis were similar to that observed for functional divergence; however, the potential peptide-binding repertoire was very similar for all breeds. These results provide an insight into MHC evolution and contribute toward efforts at describing bovine MHC diversity according to breed and location.

Study Population and DNA Extraction
Whole blood was collected aseptically from the coccygeal vein of 165 Normande and Morucha cattle. Considering that Colombian and Spanish cattle farming is characterized by extensive livestock production systems having a relatively few animals per herd, phenotypically purebred animals were sampled from the regions having the greatest amount of purebred cattle in Colombia (Normande) and Spain (Morucha), sampling the most representative amount of farms from each region ( Table 1). The herds and purebred animals analyzed were sampled randomly, avoiding related individuals/farms for each region. Normande cattle samples came from Colombia's Nariño (n = 30), Cundinamarca (n = 41), and Boyacá departments (n = 40). Morucha cattle samples (n = 54) came from Spain's Salamanca province. A PureLink Genomic DNA Mini Kit (Invitrogen, Carlsbad, CA, USA) was used for extracting genomic DNA (gDNA), following the manufacturer's instructions. Previously reported data regarding further 17 taurine and zebu populations from Asia and South America was also included for comparison (Takeshima et al., 2003;Giovambattista et al., 2013;Takeshima et al., 2015a;Takeshima et al., 2015b;Takeshima et al., 2018). Table 1 gives general information about the populations analyzed. The Universidad de Ciencias Aplicadas y Ambientales (UDCA, Bogotá) Research Ethics Committee (Minute No. 201901) approved this study.

Polymerase Chain Reaction, Cloning and Sequencing
Sequence-based typing was used for identifying BoLA-DRB3 alleles. The BoLA-DRB3 gene's exon 2 was amplified by polymerase chain reaction (PCR), using a high-fidelity polymerase along with DRB3F and DRB3R primers (Ledwidge et al., 2001). The PCR reaction mixture contained 1X Pfx amplification buffer, 300 mM of each dNTP, 0.45 mM of each primer, 1 mM MgSO 4 , 1 U Platinum Pfx DNA Polymerase (Invitrogen), and 50 ng gDNA in a 50 ml final volume. Two independent reactions were performed for each sample, following Lenz and Becker's recommendations (Lenz and Becker, 2008) to avoid the formation of chimeric products; the thermal profile was as follows: a denaturation step at 94°C for 5 min followed by 30 cycles of 94°C for 30 s, 64°C for 30 s, and 68°C for 1 min. The Wizard SV Gel and PCR Clean-Up System (Promega, Madison, WI, USA) was used for purifying PCR products which were then sequenced using a BigDye Terminator kit (Macrogen, Seoul, South Korea and Sequencing Service of the University of Salamanca, Spain) in both directions (using the DRB3F and DRB3R primers).
Another two PCR products were purified and ligated independently into pELMO vector for those animals to which  (Takeshima et al., 2003) a genotype could not be assigned from the PCR products sequences (Ramos et al., 2017). The ligated PCR products were used for transforming Escherichia coli TOP10 chemically competent cells (Invitrogen); colony PCR, with the DRB3F and DRB3R primers, was used for screening the resulting colonies for positive recombinant clones. The colony PCR reaction mixture contained 1X NH 4 reaction buffer, 250 mM of each dNTP, 0.25 mM of each primer, 1.5 mM MgCl 2 , 1 U BIOLASE DNA polymerase (Bioline, London, UK) in a final 10 ml volume. A Zyppy Plasmid Miniprep Kit (Zymo Research Corporation, Irvine, CA, USA) was used for plasmid extraction, following the manufacturer's instructions. A BigDye Terminator kit (Macrogen) was used for sequencing clones in both directions, using ccdBsec-F and ccdBsec-R external primers (Ramos et al., 2017). At least eight clones were sequenced for each animal.

Sequence Analysis
CLC Main Workbench v.3.6.5 software (CLC bio, Aarhus, Denmark) was used for assembling each independent PCR and clone sequences (manually edited if necessary) and heterozygous positions were identified for producing a consensus containing IUPAC ambiguity codes. The genotype for each animal was assigned by comparing these sequences with BoLA-DRB3 allele sequences reported in the IPD-MHC database (Maccari et al., 2017) using HAPLOFINDER (Miltiadou et al., 2003), according to Baxter et al., 2008.

Genetic Diversity Measurements, Hardy-Weinberg Equilibrium and Selection
Allele frequencies and the amount of alleles (N a ) were obtained by direct counting; 95% confidence intervals (95%CI) were calculated for allele frequencies, according to Fung and Keenan (Fung and Keenan, 2014), allowing for a small sample and population size. Observed heterozygosity (h o ) and unbiased expected heterozygosity (h e ) according to Nei (Nei and Chesser, 1983) were estimated using Arlequin v.3.5 software for population genetic analysis (Excoffier and Lischer, 2010). Five randomized datasets of 30 animals were used to calculate h o to confirm that no sampling size effect would alter diversity estimates. F IS statistic (the correlation of alleles within an individual relative to the subpopulation) (Weir and Cockerham, 1984) included in Genepop v.4.7.0 was used for estimating potential departures from Hardy-Weinberg equilibrium, using the exact test of significance (Rousset, 2008). Arlequin v.3.5 was used for calculating the amount of polymorphic sites (S), nucleotide diversity (p), and the average amount of pairwise nucleotide differences between populations (Nei and Li, 1979). MEGAX software (Kumar et al., 2018) was used for calculating the average amount of synonymous (d S ) and nonsynonymous (d N ) substitutions per site by Nei-Gojobori's method with Jukes-Cantor correction and Z-test was used for assessing d N /d S ratio significance (Nei and Gojobori, 1986).

Population Structure and Differentiation
Normande population structure was evaluated first because samples from this breed were obtained from three geographically defined regions having local breeding practices which could have differentiated cattle populations. Population structure and genetic differentiation between populations were evaluated by estimating F ST statistics (the correlation of randomly chosen alleles within the same subpopulation relative to the entire population), as described by Weir and Cockerham (Weir and Cockerham, 1984) using Arlequin v.3.5 (Pairwise F ST ) and Genepop v.4.7.0. (overall F ST ). POPTREE2 (Takezaki et al., 2010) was used for estimating genetic distances D A (Nei, 1978) from allele frequencies and constructing dendrograms using the NJ algorithm (Takezaki et al., 2010). PAST software v.3.2 (Hammer et al., 2001) was used for assigning confidence levels to branch nodes by bootstrapping D A genetic distances, with 10,000 replicates. PAST software v.3.2 was used for metric multidimensional scaling (MDS) analysis, based on D A genetic distances, and principal component analysis (PCA), based on allele frequencies.

PBR Similarity and Covariation
GeneDoc (Nicholas et al., 1997) was used for calculating identity and similarity percentages (assessed by the BLOSUM62 substitution matrix) for the 31 putative positions constituting the MHC-DRB PBR (Suarez et al., 2006). Percentage similarity was calculated for all observed genotypes within populations; five randomized datasets of 30 individuals from each population were used for one-way ANOVA and Bonferroni tests, after using STATA for proving normal distribution and homoscedasticity (StataCorp, 2017).
WebLogo software (Crooks et al., 2004) was used for creating a logo of PBR positions for each population, using the BLOSUM62 substitution matrix, including a color similarity scheme. Each position's frequency depended on the allele frequency observed for each population, considering only those having greater than 5% frequency. Covariation was estimated by constructing a vector having 620 coordinates (31 positions × 20 possible aa) for each population. Another Multidimensional Analysis Package (amap R 3.4.3) (R_Core_Team, 2019) was used for calculating Pearson correlation coefficients (PCC) between vectors. The R 3.4.3 stats package (R_Core_Team, 2019) was used to perform a metric MDS on the correlation distance matrix so obtained.

Locus and Hardy-Weinberg Equilibrium Regarding Normande and Morucha Breeds
Overall F ST (F ST = 0.022) and pairwise F ST calculated within the Normande breed indicated that the three populations were differentiated at a similar degree to that seen between breeds, and therefore these three Colombian Normande populations were considered separately for diversity analysis.
We determined the observed (h o ) and expected (h e ) heterozygosity and the fixation index (F IS ) in Morucha and Normande cattle populations and compared them to those reported for other cattle breeds to examine BoLA-DRB3 locus variability and potential departures from Hardy-Weinberg equilibrium. No relation with the original sample size was detected when h o was calculated from randomized datasets of 30 animals, indicating that this estimate (and its derivatives) can be used to compare genetic diversity in populations having unequal sized samples. For the three Normande populations as well as the Morucha population, h e was higher than h o with a significant heterozygote deficiency ( Table 3). All three Normande populations as well as the Morucha had the lowest h o compared to the other cattle populations, excluding the Chilean Hereford-Jersey population, although having similar h e levels. These populations also had the greatest deviations from Hardy-Weinberg equilibrium (F IS ).

Genetic Diversity at Sequence Level and Selection Pattern in BoLA-DRB3
The amount of polymorphic sites (S) and nucleotide diversity (p) were estimated to evaluate nucleotide genetic variability within populations ( Nucleotide genetic variability between populations was also assessed calculating the average amount of pairwise nucleotide differences ( Figure 1A; Supplementary Data Sheet 2), which ranged from -0.0367 (between Chilean Holstein and Peruvian Holstein) to 4.73 (between Chilean Hereford-Jersey and Peruvian Nellore-Brahman). Nucleotide differences between populations can be negative in rare cases where nucleotide diversity is high because it implies subtracting variability within populations from total variability (Nei and Li, 1979).

Population Structure and Genetic Differentiation
Pairwise F ST values were used for assessing genetic differentiation between populations ( Figure 1B; Supplementary Data Sheet 3). Pairwise comparisons ranged from −0.0010 (Chilean Holstein with Peruvian Holstein) to 0.1338 (Chilean Hereford with Peruvian Nellore-Brahman). Similar to that seen regarding the amount of pairwise nucleotide differences, the lowest F ST values were observed when comparing Holstein populations to one another. Nevertheless, except for Bolivian Gir, zebu populations had marked differentiation compared to taurine populations. Similar to that seen in Chilean Hereford and Japanese Jersey, Spanish Morucha was differentiated more from other taurine  populations than Normande (0.0688 average pairwise F ST in Spanish Morucha and 0.0470 in Normande). Nariño Normande was also observed to be most divergent from zebu populations and Boyacá Normande was less differentiated from taurine populations than Boyacá Normande and Nariño Normande. The highest differentiation was observed between Bolivian Nellore with Holstein populations and when comparing a group consisting of Japanese Jersey, Spanish Morucha, and Chilean Hereford to Bolivian Nellore and Peruvian Nellore-Brahman. The remaining comparisons had an intermediate level of differentiation. The D A genetic distance matrix was used to construct a dendrogram, using the NJ algorithm (Figure 2A and Supplementary Data Sheet 4). D A genetic distances were also used for depicting the relationship between populations in two dimensions by MDS analysis (Figure 2B). Dendrograms and MDS gave a similar distribution pattern, only differing regarding the clustering of groups not well supported by bootstrap (Japanese Black, Japanese Shorthorn, Boyacá Normande, Chilean Hereford-Jersey, Bo_Yac, Chilean Wagyu, Colombian Hartón del Valle and Chilean Hereford). MDS analysis had a 0.799 R-square (RSQ), indicating that representation in two coordinates gave a good description of real BoLA-DRB3 variability. Zebu populations were clearly separated from taurine populations. Despite the two types of cattle having a heterogeneous distribution, lacking any clear clustering except for the Holstein group, we were able to identify four clusters which were relatively well supported by bootstrap (>90). The zebu group (Bolivian Nellore, Peruvian Nellore-Brahman and Bolivian Gir), the Holstein group (Japanese Holstein, Peruvian Holstein, Peruvian Holstein, Chilean Holstein, Bolivian Holstein and Argentinian Holstein), a group formed by Cundinamarca Normande and Nariño Normande and a group formed by Spanish Morucha and Japanese Jersey. Boyacá Normande was noted to be closer to the other taurine populations than Cundinamarca Normande and Nariño Normande.
Allele frequencies were used for PCA and these results were compared to MDS to identify which alleles contributed most to differentiate or cluster cattle populations. Groups clearly identified in both MDS and in the first two PC were depicted as having the major alleles contributing to such grouping ( Figure  2B and Supplementary Data Sheet 5).
Population distribution in the first two principal components (PC) (Supplementary Figure 1), accounting for 44.6% of total variance, was very similar to MDS. Boyacá Normande was located toward the Holstein variability region in the first PC and was clearly differentiated from Spanish Morucha, and the other Normande, which were located toward zebu populations' variability region. Boyacá Normande differentiation in the first PC was explained by the high BoLA-DRB3

PBR Sequence Similarity and Covariation Between Populations
PBR position similarity was calculated for the alleles reported in the populations used in this study and from previous reports (Takeshima et al., 2003;Giovambattista et al., 2013;Takeshima et al., 2015a;Takeshima et al., 2015b;Takeshima et al., 2018) to further analyze BoLA-DRB3 diversity ( Five randomized sets of 30 animals were constructed to reduce possible sampling effects, considering the clustering observed in diversity analysis results. The groups compared were Holsteins, zebu populations, Nariño Normande, Cundinamarca Normande, Boyacá Normande, and Spanish Morucha. Significant differences between the evaluated groups were found for all randomized datasets (≤0.05 p-value). The Bonferroni test gave significant differences between Spanish Morucha and Holsteins populations (≤0.05 p-value) for all randomizations. Holstein populations' percentage similarity was lower than for Spanish Morucha (81.61%), indicating that Spanish Morucha had less diversity in the PBR. When identity in PBR positions was assessed, we found that 129 out of 140 alleles had a different PBR.
The logo representation of PBR had remarkable similarities among populations, having very similar substitution patterns ( Figure 3A). PBR logos had highly variable positions (11, 13, 37, 70, 71, and 74), some of which were shared by several pockets and usually had nonconservative substitutions. Covariation analysis of PBR aa sequence logos for all cattle populations ( Figure 3B) showed that most variability was due to differences in aa frequencies and not aa variations in the same position. Global covariation between PBR positions was high (PCC = 0.92). Three groups could be discriminated (Figures 3B, C). The first contained Cundinamarca Normande, Nariño Normande, Chilean Hereford, Japanese Jersey, Chilean Hereford-Jersey, and Spanish Morucha (group PCC = 0.94). The second included Japanese Shorthorn, Chilean Wagyu, Japanese Black, Boyacá Normande, Paraguayan Holstein, Argentinian Holstein, Japanese Holstein, Peruvian Holstein, Chilean Holstein, Bolivian Holstein, Bo_Yac, Bolivian Gir, and Colombian Hartón del Valle (group PCC = 0.96) and the third consisted of Peruvian Nellore-Brahman and Bolivian Nellore (group PCC = 0.93). All taurine populations and zebu Bo-Gir had an average 0.94 PCC. Although these groups were different regarding previous analysis, similar groups were found, distinguishing zebu populations (Peruvian Nellore-Brahman -Bolivian Nellore), Holstein populations and the Spanish Morucha -Japanese Jersey and Cundinamarca Normande -Nariño Normande association. Overall covariation gave Spanish Morucha as the most differentiated population between taurine breeds, similar to zebu Bolivian Nellore and Peruvian Nellore-Brahman.

DISCUSSION
MHC gene analysis can provide information about domestic species' level of overall genetic diversity. While several studies have shown decreased MHC variability in association with population bottlenecks (Bollmer et al., 2007;Bollmer et al., 2011;Mason et al., 2011;Zhang et al., 2016), others have found that high variability persists despite extreme population bottlenecks (Edwards and Hedrick, 1998;Garrigan and Hedrick, 2001;Aguilar et al., 2004;Borg et al., 2011;Moutou et al., 2013;Newhouse and Balakrishnan, 2015). The BoLA-DRB3 gene is the most polymorphic bovine MHC class II DR gene and is the only one having been shown to be fully functional, although the total amount of alleles is low compared to species such as humans (Shiina et al., 2017). This study was therefore aimed at describing BoLA-DRB3 genetic diversity in two new cattle breeds and analyzed whether domestic cattle breeding has decreased actual MHC variability.
Several measures of genetic diversity, such as h s or F IS , have corrections for sample size (Nei and Chesser, 1983). As expected, h o estimated for random datasets did not depend on original sample size. This was because h o is an unbiased estimator of parametric value (Nei and Chesser, 1983) and is mainly determined by the sampling method used. Therefore, the lower h o and higher F IS value for Spanish Morucha and Colombian Normande compared to the other cattle populations (Takeshima et al., 2003;Giovambattista et al., 2013;Takeshima et al., 2015a;Takeshima et al., 2015b;Takeshima et al., 2018) indicated heterozygote deficiency for these two breeds.
Drift and selection interact to configure the allele distribution observed in populations. MHC (a nonneutral marker) allele distribution would have been expected to be that of balancing selection with excess heterozygotes (Hughes and Nei, 1988;Hughes and Yeager, 1998;Takeshima et al., 2008). Nevertheless, the heterozygote deficiency observed in Spanish Morucha and Colombian Normande suggested that other evolutionary processes, such as inbreeding and/or bottlenecks, would have been acting on these populations. Genetic improvement regarding dairy (Normande) or beef production (Morucha) traits or bottlenecks in these breeds' origin might have contributed to the low diversity observed and high level of inbreeding inferred (>0.2 in Normande and Morucha) for the BoLA-DRB3 locus. Decreased MHC diversity can result in a narrower spectrum of pathogens which a population can recognize, less viability and even loss of pregnancies due to maternal-foetal interactions (Sommer, 2005). However, functional diversity could be preserved despite few alleles being involved.
BoLA-DRB3 divergence between Colombian Normande populations was even higher in some cases than the divergence observed between breeds on different continents (e.g., Spanish Morucha and Japanese Jersey). BoLA-DRB3 divergence between Normande populations can be highlighted compared to that for Holstein populations which are not so differentiated, despite also being phenotypically homogeneous. Comparing French and Colombian Normande cattle might provide insights into Normande populations' regional divergence because the French Normande breed is slightly variable and is considered a completely closed breed (Danchin-Burge et al., 2012;Porter et al., 2016).
Spanish Morucha is an autochthonous cattle breed from the black Iberian group which is restricted to the Northwest of the Iberian Peninsula; it was initially used for draft, milk and beef but is now solely used for meat production. It has moderate inbreeding rates and a population structure restricting mating within herds (Canas-Alvarez et al., 2014). Nevertheless, a restricted gene flow regarding other autochthonous Spanish breeds having geographical proximity during their history and similar production systems might have moulded Spanish Morucha genetic diversity (Canas-Alvarez et al., 2015). As such, this breed has had relatively low human intervention and can thus be considered close to a natural cattle population.
Spanish Morucha has been shown to be more differentiated from other taurine and zebu breeds. Likewise, Japanese Jersey and Chilean Hereford pairwise F ST values for Spanish Morucha were the highest among taurine populations. A similar trend was found when analyzing the sequence information for each allele through the average amount of pairwise nucleotide differences, Spanish Morucha and Japanese Jersey being more differentiated. BoLA-DRB3*048:02 high frequency was the main reason for Spanish Morucha and Japanese Jersey clustering. As European cattle were introduced to improve Japanese native breeds, a similar BoLA-DRB3 gene pool in Iberian and Jersey founder populations or the direct introduction of Iberian cattle along with similar selective pressures might have led to these two cattle populations' characteristic divergence, despite these two breeds having different production features and purposes.
The frequency of an MHC allele or group of alleles might become increased in cattle populations having the same selective pressures imposed by similar production systems, even in different geographic locations. From the analysis of previously published data (Takeshima et al., 2015a), we found that Holstein populations were enriched with alleles reported to be associated with resistance (BoLA-DRB3*011:01 and BoLA-DRB3*012:01) and susceptibility (BoLA-DRB3*001:01 and BoLA-DRB3*015:01) to mastitis (Yoshida et al., 2009b;Yoshida et al., 2009a). The Holstein breed has been submitted to intense selective pressure due to milk production traits which can make it susceptible to mastitis but, at the same time, there is also pressure for selecting resistant individuals. These two forces might thus create high frequencies of alleles associated with susceptibility and resistance to be in balance. However, it is worth noting that Japanese Black (a beef breed) was also enriched by these alleles. No particular enrichment with this kind of alleles was observed for Normande or Morucha cattle.
The greatest variability in the b1 domain is located in PBR positions, thereby determining the peptides to which an allele can bind. In agreement with other diversity analysis in which Spanish Morucha had lower nucleotide and codon variability, variability in the PBR was lower and could indicate that the magnitude of this population's peptide repertoire could be significantly smaller than that of other cattle. Holstein populations shared similar PBR constitution and covariation patterns by contrast with Normande ones, which had a more differentiated PBR. Despite zebu populations being welldifferentiated, Bolivian Gir appears to have a BoLA-DRB3 allele distribution similar to that of taurine breeds, thus grouping them with taurine populations. High correlation coefficients in PBR aa sequence logo analysis indicated that all breeds tended to have similar peptide repertoires, despite genetic differentiation. This suggested that domestication and/or breeding have not decreased functional MHC variability. Indeed, the fact that Spanish Morucha has lesser variability suggests that breeding might have increased both genetic and functional variation.
Great diversity was observed when analysing BoLA-DRB3 identity in the PBR, despite having fewer alleles compared to other species, such as primates (Suarez et al., 2006), i.e., even though having fewer alleles, potential peptide repertoire size in bovines is similar to that for other species. This means that cattle have 129 distinct PBR (based on sequence identity) and that infectious disease control programmes should be designed as if there were 129 different alleles instead of 140 (136 in the IPD Database and four reported elsewhere). A similar amount of PBR sequences has been determined in humans and owl monkeys (Suarez et al., 2017). This is especially important when it comes to peptide-based vaccine design which might be more efficacious when the peptide-MHC complex has potentially greater affinity as this can be increased by modifying peptides to bind to MHC molecules [reviewed in (Patarroyo et al., 2011)].
However, allele frequency would be the main guide regarding population control programmes based on MHC diversity. For example, based on allele frequency distribution, 67% to 73% of the peptide repertoire in Holstein populations would be determined only by seven alleles and this information could be used when designing peptide-based vaccines. A potential limitation of such approach is that the MHC-DR PBR positions used have been documented by X-ray crystallography just in humans and mice. Nevertheless, several of these positions have been associated with vaccine responses and susceptibility or resistance to infectious diseases (Garcia-Briones et al., 2000;Baxter et al., 2009;Rastislav and Mangesh, 2012;Carignano et al., 2017) and thus analysis based on these positions is reasonable.

CONCLUSIONS
Despite genetic divergence being observed between cattle populations, their peptide repertoires seem to be functionally equivalent. The great similarities and covariation among populations regarding the BoLA-DRB3 PBR might signify that BoLA-DRB3 peptide repertoire size is very similar among the cattle populations analyzed here and thus breeding has not reduced cattle's MHC functional diversity. Although cattle have fewer alleles than species such as some primates, PBR sequence identity supports the idea that the potential peptide repertoire between these species is equivalent and functional diversity can persist despite population bottlenecks. In other words, less genetic MHC diversity (fewer alleles) does not mean loss of functional diversity. Genetic improvement regarding breeding and domestication appears to have produced three different Colombian Normande populations based on BoLA-DRB3 genetic diversity. When analysing these populations regarding the diversity described for other cattle, it was noted that the forces which have contributed toward moulding BoLA-DRB3 locus diversity in Colombian Normande have produced remarkable differentiation compared to other phenotypically homogeneous breeds, such as Holstein. It is also remarkable that Spanish Morucha was shown to be strongly differentiated in all types of analysis.

DATA AVAILABILITY STATEMENT
All datasets generated for this study are included in the article/ Supplementary Material.

ETHICS STATEMENT
The animal study was reviewed and approved by The Universidad de Ciencias Aplicadas y Ambientales (UDCA, Bogotá) Research Ethics Committee (Minute No. 201901) approved this study.

AUTHOR CONTRIBUTIONS
MB, DO, CS, and MP designed the study. MB, DO, and CV carried out experimental procedures. MB, DO, CS, BV, CV, JL-A, AM, IO, and MP analyzed the data and wrote the manuscript. AM and MP supervised the study. All authors approved the final version of the manuscript.

FUNDING
This study was partially supported by the Diputación de Salamanca, Caja Rural de Salamanca and the Instituto de Salud Carlos III, ISCIII, Spain (www.isciii.es) under grants: RICET RD16/0027/0018 and PI16/01784. MB was partially financed trough a scholarship from the Microbiology Postgraduate Program, Universidad Nacional de Colombia. DO was financed through a scholarship from the PhD Program in Biomedical and Biological Sciences, Universidad del Rosario.

ACKNOWLEDGMENTS
We would like to thank Doctor Alejandro Garavito and the Colombian Normande Breeders Association for their support in obtaining Cundinamarca Normande samples and José Manuel Sánchez Recio from the Genealogic National Morucha Breeders Association for obtaining the Spanish Morucha samples. We would also like to thank to Jason Garry for extensively reviewing the manuscript.