BLG: β-Lactoglobulin; bp: Base-Pair; LD: Linkage Disequilibrium; R: Restriction Site for RsaI ; ‘+1’: Transcription Start Site; TRE: Thyroid Receptor Element; PRE-RC: Progesterone Receptor Element-Reverse Compliment; AP: Activator Protein; NF1: Nuclear Factor 1; MAF: Mammary Cell Activating Factor; MGF: Mammary Gland Activating Factor; MPBF: Milk Protein Binding Factor; PCR-RFLP: Polymerase Chain Reaction-Restriction Fragment Length Polymorphism

Bovine milk contains 3-5% protein, of which 80% is casein and 20% is whey protein [1]. These proteins play a vital role in coagulation and curdling of milk. Among the whey proteins, β-lactoglobulin (BLG ) forms the major portion in the milk of ruminants. This gene has been mapped to bovine chromosome 11 [2]. Studies of BLG gene in taurine cattle have revealed 11 alleles [3] and the most frequent alleles A & B being are observed to be associated with milk composition and cheese making properties [4]. Allele A of BLG encodes amino acids aspartic acid and valine, while allele B encodes Glycine and Alanine at positions 64 and 118, respectively [3,5]. In bovine milk, BLG protein variants A and B are associated with different amounts of BLG protein. The variant A has a higher BLG protein concentration than variant B [6,7]. In past, several groups have reported the significant effect of protein variants A and B on milk protein composition or milk production traits [8-15]. These studies revealed that BB genotype is associated with higher casein and lower BLG content. Comparatively, the B variant of BLG was observed to be expressed at relatively lower level in milk with a concomitant increase in casein number [16-18]. In a comprehensive study encompassing 8000 animals and three lactation records, Ng-Kwai-Hang et al. [9] observed that with replacement of allele A by B, there was a significant increase in protein percentage in all the three lactations whereas protein yield increased significantly only in the second lactation [9].

In another study involving 6803 animals with first lactation records, no significant effect on protein percentage was observed in animals with BB genotype [19]. However, the effect of BB genotype on milk yield (p=0.019), fat% (p<0.001) and protein yield (p=0.009) was significant. Sitkowska et al. [20] reported higher milk yield, higher fat and protein content in animals with AA genotype, however the study was based on only 155 cow and also significance of the effect was not provided [20]. Ozdemir and Dogru [21] found significant effect of AA genotype on increased fat percentage and AB genotype on increased milk yield [21]; while Kucerova et al. [22] observed non-significant effect of AA/AB/BB genotypes on milk traits [22].

Efforts have also been made to analyze the variations in the BLG promoter and annotate transcription factor binding sites for milk protein gene in taurine cattle [8,23-27]. Wagner et al. [23] characterized the promoter region of Holstein Friesian cattle and identified a total of 16 alleles (R1A/B to R16A/B) [23]. Amongst these, allele R5A/B at position -204 (G>C) from the transcription start site holds importance as it lies within the binding site of milk protein binding factor (MPBF), and might affect the activity of the gene product. The polymorphism analysis of bovine BLG promoter region by Lum et al. [24] revealed 10 polymorphic sites [24]. They confirmed functional importance of transversion (G to C) within a consensus binding site for activator protein-2 (AP-2) at position –430 bp from the transcription initiation site. Wagner et al. [21] also confirmed the association of alleles (A/B) in the coding region with variations in the promoter region [23]. Such associations might explain the difference of protein variant dependent BLG synthesis observed in vivo .

In contrast to taurine cattle, no such reports are available on polymorphism analysis of BLG gene in diverse Indian cattle breeds. Apart from few studies [28,29] wherein protein variants A and B were genotyped using PCR-RFLP in a small number of breeds, no other studies have been carried out to analyze polymorphism in the coding and promoter region. The present study was therefore planned with the objective to analyze distribution pattern of allelic variants in the coding (protein variants) and promoter region across diverse Indian native cattle breeds.

Cattle breeds and extraction of genomic DNA

A total of 968 animals belonging to 23 Indian native cattle breeds were genotyped for BLG protein variants A, B and R5 or -204 G/C of BLG promoter region. Out of these, 779 animals were analyzed for the BLG promoter region. These animals were subset of the 968 animals analyzed for exon4 region. Animals that were not genotyped for promoter region were considered missing values and excluded in verification of linkage between the variants of the two analyzed regions. The details of cattle breeds, their geographical location, agro climatic zones, utility and coat color is presented in Table 1. Blood samples were collected from unrelated and true-breeding animals randomly picked from the respective breeding tracts of different cattle breeds using vacutainer tubes containing EDTA as anticoagulant. Genomic DNA was isolated using standard proteinase K digestion followed by phenol: chloroform extraction procedure [30].

Table 1: List of Indian native cattle breeds included in the study and their categorisation based on utility type, geographical location, coat color and agroclimatic zone along with milk fat%

PCR-RFLP genotyping of BLG protein variants

A total 968 animals representing 23 cattle breeds were genotyped to uncover the status of one of the SNPs associated with BLG protein variants A and B. Frequency distribution of A and B alleles was analyzed within a 247 bp fragment, covering exon 4 and flanking intron 4 region of BLG gene with restriction site for HaeIII , using the primer pairs reported by Shetty et al. [31]. For variant A, the restriction site for HaeIII exists within intron 4 only, whereas, for variant B restriction site exists both in exon 4 and intron 4. The PCR reaction was carried out in a reaction volume of 25 μl containing 100 ng of genomic DNA, 5 pmol of each primer, 1.5 mM MgCl2, 200 μM dNTP and 1.0 unit of Taq polymerase (MBI Fermentas). Thermal cycling condition for amplification of the coding region was 95°C for 5 min, followed by 30 cycles at 94°C for 30 sec, 59°C for 40 sec and 72°C for 20 sec with a final extension at 72°C for 3 min. For the excision reaction, 10 μl of each PCR amplified product was digested overnight with 3 units of HaeIII restriction endonuclease (MBI Fermentas, USA) at 37°C. Restriction fragments were separated by electrophoresis on ethidium bromide stained 3% agarose gel in 1x TAE buffer and visualized using gel documentation system (UVP, Cambridge, UK).

Sequence and PCR-RFLP analyses of BLG promoter region

A 705 bp fragment comprising of 584 bp of the proximal promoter region and 121 bp of exon 1 was amplified using the oligonucleotide primers described by Yahyaoui et al. [32]. The composition of PCR reaction mixture was as described for the coding region. A touchdown PCR protocol was employed with temperature profile of 95°C for 5 min, 10 cycles of 97°C for 15 sec, 63°C for 1 min, and 72°C for 90 sec followed by 25 cycles of 95°C for 30 sec, 63°C for 1 min and 72°C for 90 sec followed by a final extension at 72°C for 5 min. Purified DNA sample was bidirectionally sequenced with ABI Prism® BigDyeTM Terminator Cycle Sequencing Kit (version 3.1) (Applied Biosystems, Foster City, CA).

^aNR, northern; NWR, northwestern; CR, central; SR, southern; ^bR, red; G, gray; W, white; B, brown; ^cAgroclimatic zones: SAR, semiarid; AR, arid; TWDR, tropical wet and dry; HSTR, humid subtropical; ^d(Nivsarkar et al., [41]); ^ehttp://dad-training.fao.org; NA: Not available.

All sequencing reactions were analyzed using an ABI 3100 DNA capillary sequencer (Applied Biosystems, Foster, CA). Sequence comparisons were done using Phrap and Gap4 integration from the Staden Package, a suite of sequence handling and analysis software developed at the Medical Research Council Laboratory of Molecular Biology (Cambridge, UK). The sequences were annotated for binding site of transcription factors using program MATCH v1.0, AliBaba2.1, TESS and information available in literature [23,24,27]. The identified SNP (R5A/B, -204 G/C) was further examined using PCR-RFLP on all the 779 animals (subset of 968 animals genotyped for protein variant A and B) analyzed for the promoter region. Restriction digestion was carried out using 10 μl of each PCR product and 3 units of RsaI enzyme (MBI Fermentas, USA) overnight at 37°C. The digested products were resolved on 3% agarose gel in 1x TAE buffer and visualized under UV trans illuminator

Statistical analysis

Gene and genotypic frequency was calculated by direct counting method. Chi-square test was carried out manually to test the populations for Hardy-Weinberg equilibrium as well as for differences between the animals on the basis of geographical location, agro climatic zones, utility and coat color. Linkage disequilibrium (LD) measures, D’ and r2, between all variations within promoter region and two protein variants were estimated using Arlequin v3.5 software [33].

Detection of BLG protein variants A and B

Screening of the exon 4 region of BLG using HaeIII detected protein variants A and B with three different genotypes (AA, AB and BB). The restriction fragments corresponding to AA, AB and BB genotypes were: 148 and 99 bps; 148, 99 and 74 bps, and 99 and 74 bps, respectively (Figure 1a). The gene and genotype frequencies of BLG protein variants A and B in each of the analysed Indian cattle breeds are depicted in Table 3. The genotype BB was predominant in 15 breeds and ranged from 0.214 (Malnad Gidda) to 0.877 (Kangayam), with an average frequency of 0.531. On the other hand, AB genotype was in higher proportion in 7 breeds with an average frequency of 0.406. The frequency of AB genotype ranged from 0.069 (Ponwar) to 0.755 (Mewati). The least common genotype, AA with an average frequency of 0.064, was absent in 5 breeds viz . Hariana, Mewati, Kangayam, Kherigarh and Red Sindhi.

Figure 1: Agarose gel showing BLG (a) exon 4-intron 4 and (b) promoter variants of Indian native cattle breeds; M: 100 bp DNA ladder. PCR-RFLP of the coding region detects one of the SNPs associated with BLG protein variants A and B, while the promoter region detects the R5 or -209 G/C variants.

Each of the studied cattle breeds demonstrated predominance of allele B with an average frequency of 0.740. The frequency of allele B ranged from 0.559 (Gir) to 0.938 (Kangayam). On the other hand, average frequency of allele A was 0.260, with Gir showing the highest frequency (0.441). The relatively high frequency of B allele in Kangayam corresponded to higher prevalence of BB genotype in the breed.

The presence of allele B at high frequency and predominance of BB and AB genotypes across the Indian native cattle breeds is similar to earlier reports on Indian cattle. Ganai and Bhat [34] reported high frequency of allele B (82.9%) and BB genotype (71.4%) in Hariana cattle. Similarly, Rachagani et al. [33] reported high prevalence of B allele (0.83 and 0.61) and BB genotype (0.693 and 0.244) in Sahiwal and Tharparkar cattle, respectively. In taurine cattle, comparatively lower frequency of B allele and BB genotype was reported in Holstein (0.647 and 0.399 [11]; 0.471 and 0.198 [34]), Czech Fleckvieh (0.489 and 0.240), Brown Swiss (0.56 and 0.296) [21] and Jersey (0.66 and 0.392) [31]. The allelic distribution revealed the higher frequency A allele in taurine breeds and conversely the naturally evolved Indian zebuine cattle breeds have higher frequency of B allele. The lower occurrence of BLG protein variant B and higher occurrence of variant A allele observed in taurine cattle might be attributed due to their intensive selection for higher milk yield, as A allele and AA genotype have been associated with higher milk yield [20,35]. In contrast, the higher frequency of B allele and near absence of AA genotype breeds could be an indication for lack of any selection pressure for milk yield in the naturally evolved Indian cattle breeds.

Till now, no association data using BLG variants with dairy traits are available for Indian cattle. However, many such studies have already been reported in different taurine populations trying to evaluate the effect of BLG variants on milk production and protein content [9,19,12,20,35-37]. The effect of BLG genotypes on milk production trait as different reports indicate is not consistent. Lunden et al. [30] and Ojala et al. [31] observed no effect of BLG genotypes on milk production traits, but studies by Bovenhuis et al. [31] and Ikonen et al. [30] observed association of BLG genotypes with production traits suggesting AA genotype for higher milk production. Conversely, most of the reports consistently indicate that A variant of BLG gene is associated with a greater concentration of BLG protein [9-10,24]. In a comprehensive study on approximately 2000 animals, Heck et al. [30], confirmed association of A variant with relatively higher concentration of BLG as compared to B variant of BLG , however BLG variants did not have any effect on milk production trait. The difference in the outcome of these studies could be due to variation in the sample size or cattle types included in the study. Besides, the difference in BLG protein level of variant A and B could be due to the differential expression of respective alleles and its linkage with variations within the promoter. Defining the allelic status of A and B protein variants of BLG across Indian cattle breeds will be crucial to initiate genotype and phenotype association studies.

Polymorphism at BLG promoter region

The amplicon of the 5'- flanking region of BLG gene was sequenced and compared with Bos taurus reference sequence [NCBI: DQ489319]. The observed length of the amplicon (705 bp) comprising of 584 bp of promoter region and 121 bp of exon 1 was five base-pair shorter than the observed for goat (710bp) by Yahyaoui et al. [32] as BLG promoter in cattle has deletion of six base-pair (between +70 bp and +71bp) and insertion of single base-pair (between -227 bp and -228bp) as compared to goat sequence.

The binding sites for promoter transcription factor of milk protein genes viz . thyroid receptor element (TRE), progesterone receptor element (PRE-RC), activator protein-2 (AP-2), nuclear factor 1(NF1), mammary cell activating factor (MAF), milk protein-binding factor (MPBF), c-Rel, Elk-1 and mammary gland activating factor (MGF) of Bos indicus were similar to Bos taurus (Figure 2). Alignment of the sequences with Bos taurus reference sequence (DQ489319) revealed 15 variations at positions +117 (G/A), +99 (T/C), +78 (T/C), +36 (A/G; R1), +32 (C/T; R2), -83 (G/A), -98 (G/C; R4), -170 (A/C), -204 (G/C; R5), -216 (C/T; R6), -287 (G/A), -362 (T/C; R7), -430 (G/C; R10), -457 (C/G; R11) and -573 (A/G). The nomenclature for different SNPs (R1- R11) in the parenthesis is as suggested by Wagner et al. [31].

Figure 2: DNA sequence and annotation of the BLG promoter in Indian native cattle breeds. Variations among Indian cattle in reference to Bos taurus (DQ489319) are marked with bold and italicized nucleotide. Sequence for primer pair are indicated by underlined oligomers at 5’- and 3’-terminals

Amongst these SNPs, six variations located within the putative transcription factor binding sites (TFBSs) are +36 (A/G; AP4), -170 (A/C; c-Rel and Elk-1), -204 (G/C; MPBF), -362 (T/C; NF-1), -430 (G/C; AP-2) and -457 (C/G; TRE and PRE-RC) possibly affects binding affinity of TFs and thus milk properties [24]. The variations R1-R2, R4-R7, R10-11 observed in the present study have also been reported by Wagner et al. [31] and Ganai et al. [32] (Table 2). Reports are also available for presence of variations at positions +99 (T/C) [27], +78 (T/C) [24,26], -170 (A/C) located within the binding sites for transcription factors (TFs) c-Rel and Elk-1 [27] in taurine breeds. However, variations at positions -424 (T/G; R9) [23], -422 (G/A; R8) [23,24,26] and -22 (G/A; R3) [23] were not observed in the present study.

Table 2: Gene and genotype frequencies of BLG/HaeIII locus (BLG protein variants A and B) in Indian native cattle breeds; a: one degree of freedom; NS: non- significant; ^*: p<0.05; ^**: p<0.01

The other three variations -83 (G/A), -287 (G/A) and -573 (A/G) found in the BLG promoter region of Indian cattle have not been described in any other study. Additionally, variation at position +117 (G/A) in the exon-1 has also not been reported previously.

Amongst the variations detected in present study, change at position -204 (G/C; R5) creates restriction site for RsaI endonuclease enzyme. The frequency of this particular variation was further analysed in 779 Indian native cattle as it is located within the binding site of milk protein binding factor (MPBF), known as signal transducer and activator of transcription-5 (STAT-5). STAT-5 is a transcription activator for BLG gene in mammary cells and bind sitespecifically in the 5'-flanking region of the gene. It might be postulated that nucleotide change at STAT5 binding site could influence the activity of its gene product, which in turn could result in allele specific differences in expression of BLG gene. Similar differential expression of alleles at BLG locus has earlier been proposed to be due to SNPs in the 5'-flanking region (G/C; R5) [23,38-41].

Rsa1 digestion of the BLG promoter amplicon, involving R5 or -204 G/C polymorphism resulted in three genotypes: GG (382 and 323 bps), GC (705, 382 and 323 bps) and CC (705bp), respectively (Figure 1b). The genotype and allele frequencies of RsaI variant in BLG promoter region across Indian zebu cattle breeds are listed in Table 3. The frequency of predominant allele G ranged from 0.333 (Ongole) to 0.833 (Malnad Gidda), with a mean value of 0.639. The frequency of C allele was lower in all the breeds, except Tharparkar (0.531), Ongole (0.667) and Ponwar (0.568). Among the genotypes, heterozygous GC was observed in all the cattle breeds analyzed with average frequency of 0.444 and ranged from 0.094 (Nimari) to 0.773 (Ponwar). On the other hand GG genotype was predominant in 13 breeds and ranged from 0.000 (Ongole) to 0.733 (Malnad Gidda) with an average frequency of 0.417.

Table 3: Gene and genotype frequencies of BLG promoter/RsaI (R5 or -209 G/C) in Indian native cattle breeds;b: one degree of freedom; NS: non- significant; ^*: p<0.05; ^**: p<0.01

The genotype frequencies of BLG protein variants A and B and BLG promoter (R5 or -204) G/C polymorphism were also analyzed to test deviation from the expected Hardy-Weinberg equilibrium using Chisquare test. Results revealed presence of genetic equilibrium in majority of the analysed breeds (Table 3). However, 5 breeds (Nagori, Mewati, Ponwar, Red Sindhi and Malnad Gidda) showed deviation from genetic equilibrium for coding region and 6 breeds (Red Kandhari, Kankrej, Gir, Ongole, Nimari and Ponwar) for promoter region. The disequilibrium may reflect events such as selection, genetic drift or population subdivision as animals in these breeds were from either organized farm (Red Sindhi), low population size (Nagori, Red Sindhi) or restricted locations (Malnad Gidda, Ponwar, Mewati).

Further, in order to understand the allelic frequency distribution across different groups of cattle, breeds were classified in different categories on the basis of their utility (milch, dual and draft), geographical location (northern/north-western, central and southern region), coat color (red and grey) and agroclimatic zone (semiarid, arid, tropical wet and dry and humid subtropical) (Figure 3). Results indicated non-significant difference between the groups in both coding (BLG protein variants A and B) and promoter region (R5 or -204 G/C) as confirmed by Chi-square test (data not shown).

Figure 3: Mean frequencies of (a) BLG protein variants A and B and (b) R5 or -209 G/C variants of BLG promoter region across different groups of Indian native cattle.

Linkage disequilibrium between BLG protein variants and variation (R5 or -204 G/C) in the promoter region was also explored. However, no LD was observed between R5 and allele A/B of the coding region. Contrary to the present finding, complete LD between the BLG protein variants A and B and the R5 SNP at position -204 (corresponding to position -209 by Wagner et al. 1994 was reported by Lum et al. [24] in Holstein, Brown Swiss and Jersey cattle, and by Ganai et al. [21] in Holstein Friesian cattle. In addition, Lum et al. [24] found complete LD between the BLG protein variants A and B and SNPs R1, R2, R3, R4, R5, R6, R7, R10, R11, R13 and R14 in the promoter region.

Phylogenetic relationship for BLG promoter region among different mammalian species

Phylogenetic analysis based on UPGMA with 1,000 replicates for flanking region of beta-lactoglobulin revealed two major groups. As expectedly, Bos indicus , B. taurus , B. grunniens and B. bubalis were grouped together, while, goat and sheep formed another group (Figure 4). Analysis of genetic distance using p-distance method based on pairwise deletion revealed maximum genetic relatedness of Indian zebu cattle with Bos taurus (98.96%) followed by B. grunniens (98.27%), B. bubalis (96.20%), C. hircus (92.06%) and O. aries (91.88%).

Figure 4: Phylogenetic relationship of beta-lactoglobulin promoter region among different mammalian species.

Utilizing 968 animals representing 23 Indian native cattle breeds, the present study portrays comprehensive frequency pattern of the BLG protein variants A and B across naturally evolved Indian native cattle, wherein, systematic information of known variants at this important candidate gene was lacking. In future, correlation analysis of protein variants and variants in promoter region may be carried out to undertake genotype: phenotype association studies for dairy traits across divergent Indian cattle breeds. Further, considering the fact that B variant is predominantly present across dairy, dual and draft cattle breeds and there is no LD between the protein variants and variant R5 (-204 G/C), we might need to look out for new variation in BLG locus and or linked gene/QTL as a potential marker with regulatory effect on dairy traits in Indian cattle.

The financial support received from the Department of Biotechnology (DBT), Govt. of India, New Delhi for this work is duly acknowledged. The authors are thankful to Director, NBAGR for providing necessary facilities to carry out this work. We also thank Mrs. Pravesh K. for rendering technical support during the course of this work