Holstein cows are large animals with black and white or red and white colour pattern that produces the most milk on average per cow out of all the breeds (Fanta, 2017). Milk is very important component to prepare human food, while nursing calves depend on it as a source of nutrition (Makamu and Tyasi, 2024; Shah et al., 2024). Cow milk includes a variety of elements necessary for growth and development, including fat, lactose, proteins, vitamins and minerals (Chitra, 2022). However, these characteristics such as fat, milk yield, lactose and protein of different cattle breeds such as Jersey, Holstein Friesian, Sahiwal, Girolando cows, etc require genetic improvement (Zaglool et al., 2016; Badola et al., 2003; Barbosa et al., 2019). Most of the time, the selection of dairy cows has been based on quantitative parameters like fat, milk, or protein yield, which are thought to be regulated by several loci (Stipp et al., 2013). Therefore, genetic improvement of quantitative traits is relatively slow as productive traits can only be measured in one sex and is affected by numerous polygenes (Kusza et al., 2015). Several studies suggested that gene approach focusing on genetic variation between or within dairy cattle population is a helpful tool to find genetic markers that could be utilized in marker-assisted selection in breeding to enhance attributes related to milk production (Chessa et al., 2013; Kishore et al., 2014; Shahlla et al., 2014; Singh et al., 2015; Raschia et al., 2018). According to Ridhowi et al., (2018), it is crucial to identify markers for complex traits because marker assisted selection allows for the precise selection of particular DNA changes that have been linked to a discernible difference or impact on complex traits. Archana (2013) reported that there is a great deal of interest in using molecular genetics technology to identify certain deoxyribonucleic acid (DNA) markers associated with economically significant features to improve the effectiveness of breeding programs by selecting young buffaloes early on for use as future breeding stock. According to a study conducted by Badola et al., (2003), the beta lactoglobulin gene possesses genetic variants that are linked to Jersey, Holstein Friesian and Sahiwal cattle breeds’ milk production traits. Zaglool et al., (2016) reported that Holstein Friesian is polymorphic for beta lactoglobulin (β-LG) gene making three genotypes (AA, AB and BB) with two alleles (A and B). However, β-LG gene variation and how it relates to milk production characteristics in South African Holstein Friesian is not yet known. This study was aimed to discover the genetic markers of beta lactoglobulin gene that might be utilized as a criterion for selection during breeding for enhancement of milk production traits of South African Holstein Friesian cows. Because the application of new technologies in breeding and selection is under-utilized in South African cattle. The findings of the study could be helpful for enhancing milk production characteristics in South African Holstein Friesian cattle through marker-assisted selection and outline a strategy to bypass the time-consuming, traditional ways of selecting Holstein Friesian cattle for milk production.

Ethical approval

The norms and guidelines established by the University of Limpopo Animal Research Ethics Committee (AREC) under project number AREC/28/2023: PG were adhered to in all procedures.
 
Study area
 
The study was carried out at the Limpopo dairy farm in Louis Trichardt, South Africa’s Limpopo province. The ambient temperature around the Limpopo Dairy farm ranges between 8°C and 21°C during winter and 19°C followed by 28°C during summer (Adesoye and Dondofema, 2021). The Limpopo dairy farm is located at latitude 23°06ʹ29.1ʺS and longitude 29°50ʹ05.6ʺE.
 
Experimental and animal management
 
A total number of 100 South African Holstein Friesian cows aged 2 to 3 years provided by the farmer were used in this study. The age of the cows was determined using the farm records. The cows at Limpopo dairy farm were kept under intensive production, whereby they were kept inside a housing and given feed inside there. Clean water was also given to the cows ad-libitum. For identification of the animals, electronic ear tags were used. Automatic rotary milking parlour was the type of milking system used at the farm. Every time before milking, vetericyn utility spray and clean paper towels were used to clean the teats of the cows. The platform rotated slowly, giving the cows enough time to enter and exit at regular intervals. Once the animal entered the milking parlour there was a sensor that identified the animal and when the milking process starts, it was then recorded as per how much milk yield that animal produced that day. This information was recorded and stayed in the daily farm records. Sick and injured were not included in this study.

Milk and blood sample collection
 
Daily farm records were used to accumulate milk yield of cows. Samples of blood (5 ml) were taken underneath the tail of each cow using a butterfly needle which draws blood from the animal straight to the sterilized vacutainer tube containing EDTA as anticoagulant, which was temporarily stored inside a cooler box with gel cool packs. The blood samples were then transported to the university’s laboratory and stored there at -20°C inside a freezer until they were shipped to Inqaba Biotechnology company (525 Justice Mahommed St Muckleneuk, 200, Pretoria, South Africa). 10 ml of milk samples were taken during morning milking using a measuring cup to represent the entire milking of each cow. The samples were collected to determine the milk composition (fat%, protein%, SNF and lactose%) using an ultrasonic portable milk analyzer (milko tester model-master mini) at the Limpopo dairy farm.
 
Amplification of β-LG gene by polymerase chain reaction
 
Blood samples were sent to Inqaba Biotechnology company for the extraction of DNA and amplification by polymerase chain reaction (PCR). β-LG gene covering intron 3, exon 4 and intron 5 were amplified using a pair of forward (5ʹGCC TCA GAC TCA GTG GTGA3ʹ) and reverse (5ʹACC ACA CAG CTG GTC TCC3ʹ) primers. Primers 3.0 software and the published nucleotide sequence of the Bos taurus β-LG gene (GenBank Accession No X14710.1) were used in the creation of the primers. Twenty microliter reaction mixtures were used for the PCR reaction, each containing 10 microliters of NEB OneTaq, 2X Master mix with standard buffer, one microliter of genomic DNA (10-30 ng/μl), one microliter of each primer (10 μM) and up to seven microliters of nuclease-free water. The conditions for amplification were as follows: 5 minutes of denaturation at 94°C, 35 cycles of 94°C for 30 seconds, 50°C for 30 seconds and 68°C for 1 minute, with a final 10-minute extension at 68°C and hold at 4°C. EZ-vision® Bluelight DNA Dye-stained 1% agarose gel (Cleaver Scientific Ltd) was used to assess the molecular weight and integrity of the PCR amplicons. As follows, PCR products were cleaned using the ExoSAP protocol: 1. To prepare the Exo/SAP master mix, a 0.6 ml micro-centrifuge tube was filled with Exonuclease I (20 U/µl 50 µl) and Shrimp Alkaline Phosphatase (1 U/µl 200 µl). 2. ExoSAP Mix with 10 µl of the amplified PCR product (step 1) Mixtures in 2.5 µl were made. 3. The mixture was thoroughly combined and incubated for 15 minutes at 37°C. 4. The mixture was heated to 80°C for 15 minutes to halt the process.
 
DNA sequencing
 
PCR products of a region of introns 3 and 5 as well as exon 4 of β-LG gene were sequenced at Inqaba Biotechnology to identify the single nucleotide polymorphisms (SNPs). Following the manufacturer’s instructions, fragments were sequenced using the Nimagen BrilliantDyeTM Terminator Cycle Sequencing Kit V3.1, BRD3-100/1000 (NimaGen, Netherlands). To align the sequences, the NCBI/BLAST/blastn suit was used. After that, the ZR-96 DNA Sequencing Clean-up Kit (Zymo Research, USA) was used to clean the labelled products. Using POP7, the cleaned products were injected onto an Applied Biosystems ABI 3500XL or ABI 3730XL Genetic Analyzer equipped with a 50 cm array. The analysis of the sequence chromatogram was done with FinchTV analysis software. Sequence alignment was done using the NCBI/BLAST site.
 
Statistical analysis
 
Version 29.0 of the Statistical Package for Social Science (IBM SPSS, 2022) was utilized to analyse the data. The genotypic and allelic frequencies were computed using Hardy-Weinberg Equilibrium test software for Population Genetic Analysis. The Hardy-Weinberg theorem’s genetic equilibrium at the population level was assessed using the Chi-square test. An investigation of the marker-trait connection was conducted using the student t-test. The model that was employed was this one:

Where,
Y_ij = Phenotypic values of traits.
µ = Population mean.
G_i = Fixed effect of genotype.
e_ij = Random residual error.

Descriptive statistics
 
Descriptive summary of milk production traits for Holstein Friesian is illustrated in Table 1 below. The highest value of mean was that of milk yield per 30 days which was 997.92 L. Fat had the lowest mean value of 1.59% and the highest coefficient of variation of 32.33% whereas Solid-not-fat had the lowest coefficient of variation of 5.55%.


 
Correlation matrix
 
Phenotypical association between constituents of milk and milk yield is shown in Table 2. The findings showed that there was a highly negative significant connection (p<0.01) between milk yield per day (MYD) with protein percentage (PP) and solid-not-fat (SNF). Milk yield per 30 days (MY30D) showed highly negative statistical correlation (p<0.01) with FP and PP. There was no association (p<0.05) between MY30D and lactose percentage (LP). MY30D showed a highly significant link (p<0.01) with SNF.


 
Amplified nucleotide sequence analysis
 
Fig 1 shows amplified PCR products of β-LG gene in Holstein Friesian cows used in this study. The amplicon size of 447bp was generated during the amplification.


 
Sequence analysis and alignment on 5174T>C of exon 4
 
The analysis and alignment of sequence 5174T>C of exon 4 is shown on Fig 2. Intron 3, exon 4 and intron 5 of β-LG gene, were sequenced and two SNPs were detected on exon 4 whereas three were found on intron 3 and no SNP was found on intron 5 (Fig 2A). A polymorphism was detected with nucleotide transition from thymine (T) to cytosine (C) at position 5174 of exon 4 when compared with the β-LG gene (accession number: X14710.1). Blast was used to find the pairwise alignments of DNA as highlighted in Fig 2B. The sequence alignment results showed 5174T>C as the SNP position (red line). Blast was used to determine the protein sequence alignment as indicated in Fig 2C. The results indicated nonsynonymous SNP as highlighted with a red box. Isoleucine (I) amino acid changed to valine (V) at position 882 was found by comparing it with the β-LG gene (acc. no. Np_776354.2).


 
Sequence analysis and alignments on 5123C>G of intron 3
 
The analysis and alignment of sequence 5123C>G of intron 3 is shown on Fig 3. A polymorphism was found with nucleotide transition from cytosine (C) to guanine (G) at position 5123 of intron 3 when compared with the β-LG gene (accession number: X14710.1) (Fig 3A). Blast was used to find the pairwise alignments of DNA as indicated in Fig 3B. The sequence alignment results showed 5123C>G as the location of the SNP as highlighted in red line.



Sequence analysis and alignment on 4982G>A of intron 3
 
Sequence analysis and alignment of 4982G>A of intron 3 is shown on Fig 4. Nucleotide transition from guanine (G) to adenine (A) polymorphism was detected at position 4982 of intron 3 when compared with the β-LG gene (accession number: X14710.1) (Fig 4A). Blast was used to find the pairwise alignment of DNA as indicated in Fig 4B. The sequence alignment results showed 4982G>A as the position of the SNP as highlighted with the red line.


 
Sequence analysis and alignments on 5099T>C of intron 3
 
Gene sequence analysis and alignments on 5099T>C of intron 3 is shown on Fig 5. A polymorphism was detected with nucleotide transition from thymine (T) to cytosine (C) at position 5099 of intron 3 when compared with the β-LG gene (accession number: X14710.1) (Fig 5A). Blast was used to find the pairwise alignment of DNA. As indicated in Fig 5B. The sequence alignment results showed 5099T>C as the SNP location as highlighted in the red line.



Sequence analysis and alignments on 5251C>T of exon 4
 
The analysis and alignment of sequence 5251C>T of exon 4 is shown on Fig 6. Polymorphism was detected with nucleotide transition from cytosine (C) to thymine (T) at position 5251 of exon 4 when compared with the β-LG gene (accession number: X14710.1) (Fig 6A). The sequence alignment results showed 5251C>T as the SNP position (red line) (Fig 6B). The protein sequence alignment as determined by Blast indicated nonsynonymous SNP as highlighted in red box (Fig 6C). Furthermore, a glycine (G) amino acid change to aspartic acid (D) at 852 position was detected when comparing experimental samples with the β-LG gene (acc. no. Np_776354.2).


 
Genotypic and allelic frequencies
 
The allelic and genotypic frequencies of β-LG gene locus in the Holstein Friesian population are shown in Table 3. Two alleles and two genotypes (homozygous and heterozygous) were noted for each SNP. For SNPs 5174T>C, 5123C>G, 4982G>A, 5099T>C and 5251C>T, allelic frequencies of T, C, G, T and C were higher than that of C, G, A, C and T respectively. Genotypic frequencies of CT, CC, AG, CT and TC were higher than the genotypic frequencies of TT, GC, GG, TT and CC for SNPs 5174T>C, 5123C>G, 4982G>A, 5099T>C and 5251C>T respectively. The Chi-square (x2) test for 5123C>G showed that genotypic and allelic frequencies were not significantly different from the expectations of Hardy-Weinberg (X2 = 1.23). The results indicate a constant genotypic and allelic frequencies of population from generation to generation. However, 5174T>C, 4982G>A, 5099T>C and 5251C>T SNPs were tested and displayed incredible genetic imbalance between alleles (P>0.05). The results indicate that from generation to generation, genotypic and allelic frequencies of population changes.


 
Polymorphism information analysis
 
The genetic diversity and polymorphism information analysis of the population are shown in Table 4. The homozygosity of the beta lactoglobulin gene was higher than the heterozygosity of it for single nucleotide polymorphisms 5174T>C, 5123C>G, 4982G>A, 5099T>C and 5251C>T with effective allele number (Ne) of 1.92, 1.22, 1.92, 1.92 and 1.72 respectively. Polymorphisms information content (PIC) indicated that there were high polymorphisms within the Holstein Friesian population for SNPs 5174T>C, 4982G>A, 5099T>C and 5251C>T. However, it showed that there were moderate polymorphisms within the Holstein Friesian population for SNP 5123C>G.


 
Association analysis of β-LG gene with milk production traits
 
Marker-traits association are displayed in Table 5. The results indicated that genotypes (TT and CT), (GG and AG) and (TT and CT) of SNPs 5174T>C, 4982G>A and 5099T>C respectively, were not significantly different from MY30D, PP and SNF (p>0.05). However, they had significant difference with MYD, FP and LP (p < 0.05), with genotypes TT, GG and TT performing better than CT, AG and CT respectively for MYD. Whereas CT, AG and CT performed higher than TT, GG and TT for FP respectively.  5123C>G SNP showed significant difference between CC and CG genotypes with MY30D, FP, SNF and LP (p<0.05). Genotype CC performed well on MY30D and SNF, while CG genotype had a high performance on FP and LP. This SNP showed non-significant difference between CC and CG genotypes with MYD and PP (p>0.05). Significant difference was found between CC and CT genotypes with FP, SNF and LP for SNP 5251C>T (p<0.05), with genotype CC doing well on FP than genotype TC and TC performing better than TT for SNF and LP. However, CC and CT genotypes were not significantly different from MYD, MY30D and PP (p>0.05).


       
The ability to anticipate how one feature will change in response to selection for another makes understanding the relationships between traits vital for improving the quantity and quality of milk produced by dairy animals (El-Moghazy et al., 2015). There was a negative relationship between milk yield per day with protein % and SNF. An increase in milk yield per 30 days was noted with a decrease on fat % and protein %, but with an increase in SNF. Milk yield had no relation with lactose % and milk yield per day had no relationship with fat %. The findings of this study agree with the study that was conducted by El-Moghazy et al. (2015) who discovered that SNF was positively correlated with milk yield of Egyptian Buffaloes, however, this study also disagrees with the same study that found that fat, protein and lactose were positively correlated with milk yield. The findings of the study conducted by Alphonsus and Essien (2012) who stated that SNF, fat and protein were not significantly correlated with total milk yield of Friesian × Bunaji and Bunaji cows of Nigeria disagreed with the findings of the present study. This study agrees with the study conducted by Yoon et al., (2004) which stated that milk yield was negatively associated with protein and fat of Holstein cows in Korea. The difference between this study and other studies might be because of the use of different species, breed and environment. The result of this study implies that decreasing protein % and SNF increases milk yield per day, furthermore, fat % and lactose % does not affect milk yield per day. An increase in SNF increases milk yield per 30 days, whereas increasing fat and protein leads to a decrease in milk yield per 30 days. Lactose % does not have any effect on milk yield per 30 days.
       
Findings of the current study revealed 2 nonsynonymous novel SNPs 5174T>C and 5251C>T. This study also noted 3 other novel SNPs 5123C>G, 4982G>A and 5099T>C. A single nucleotide polymorphism (1810C>T) in exon 3 in β-LG gene of Chinese Holstein cows was discovered by Alim et al., (2015) who investigated DNA polymorphisms in the β-LG gene associated with milk production characteristics in Holstein dairy cattle in China. Mancini et al., (2013) found a SNP (C>A) at position 968 of upstream gene variant of β-LG gene on Italian Brown cattle in Italy. Yang et al., (2012) investigated polymorphism in exon 4 of β-LG gene different B precursor and its relationship with milk production traits and protein formation in Chinese Holstein and identified 3 nonsynonymous SNPs (5239C>A, 5240A>C, 5305C>T), meaning that three SNPs caused amino acid changes. Disagreement might be because of the differential expression of genes which impacts animal’s production traits. The results of this study suggest that SNPs 5174T>C and 5251C>T causes an amino acid change from isoleucine to valine and glycine to aspartic acid, respectively, which affects structure and function of the protein, meaning that the new protein formed will cause a change in the relationship between the genotypes and the traits. The population used was under Hardy-Weinberg equilibrium (HWE) for SNP 5123C>G. However, it was not under HWE for SNPs 5174T>C, 4982G>A, 5099T>C and 5251C>T. The results of the study that was conducted by Alim et al., (2015) indicated that chi-square test for SNP 1810C>T showed all genotypic frequencies in the population to fall under Hardy-Weinberg equilibrium indicating that allele frequencies stayed the same across generations. Yang et al., (2012) reported that after chi-square test the 3 SNPs (5240A>C, 5239C>A, 5305C>T) were not under Hardy-Weinberg equilibrium. This study indicate that the studied population is under HWE implying that the allelic and genotypic frequency for SNP 5123C>G of β-LG gene on Holstein Friesian cows does not change from generation to generation. However, population studied was not under HWE for SNPs 5174T>C, 4982G>A, 5099T>C and 5251C>T, implying that for these SNPs genotypic and allelic frequency changes from generation to generation.
       
Marker trait association findings for SNPs 5174T>C and 5099T>C indicated that there was no connection between genotype TT and CT with milk yield per 30 days, protein % and SNF statistically. For SNP 4982G>A genotype GG and AG had no association with milk yield per 30 days, protein % and SNF. Genotype CC and CG of SNP 5123C>G had no relationship with milk yield per day and protein %. Marker trait association results for SNP 5251C>T indicated no relationship between genotype CC and CT with protein % and milk yield. Relationship of β-LG gene polymorphism with fat, protein and milk yield in Holstein Friesian cattle in Egypt was investigated by Zaglool et al., (2016), who found 3 genotypes (AA, AB and BB) and discovered AA had more protein % and milk yield, while BB genotype recorded higher fat %, the results are not in line with the ones of this study. This study for SNPs 5123C>G on milk yield and 5251C>T on fat % agrees with the study that was done by Hristov et al., (2011) who found 2 genotypes AA and AB of β-LG gene in Bulgarian Black Pied cattle, that revealed BB genotype to have the highest effect on milk yield and fat %. For SNPs 5124T>C, 4982G>A and 5099T>C on SNF and SNP 5251C>T on milk yield, the current study agrees with that of Tolenkhomba et al., (2014) that revealed two genotypes AB and BB that had no significant impact on milk yield and SNF of Sahiwal cattle breeds of India. The results of the current study are in contradiction with the ones of the study conducted by Dogru (2015) who investigated β-lactoglobulin genetic variations in Brown-Swiss dairy cows and their relationship with quality traits and milk yield in Turkey and found no significant association between different genotypes (AA, AB and BB) of β-LG gene and milk production constituents. The difference in the current study might be due to different environmental conditions and breeds used.
       
TT genotype for SNPs 5174T>C and 5099T>C of β-LG gene might be utilised as genetic marker when enhancing milk yield per day and lactose %, whereas
genotype CT might be used to improve fat %. Genotype CC of SNPs 5123C>G might be used to increase milk yield per 30 days and SNF, while CG be used to improve lactose % and fat %. TC for SNP 5251C>T might be used as a genetic marker to increase SNF and lactose%, whereas CC be used to improve fat%. GG genotype for SNP 4982G>A of β-LG gene might be utilised as genetic marker when enhancing milk yield per day and lactose %, whereas genotype AG might be used to improve fat %.

This study conclude that protein and solid-not-fat has the potential to be used when improving milk yield of Holstein Friesian cows. Genotypic and allelic frequency for SNP 5123C>G of β-LG gene on Holstein Friesian cows does not change from generation to generation. However, allelic and genotypic frequencies change from generation to generation for SNPs 5174T>C, 4982G>A, 5099T>C and 5251C>T of β-LG gene on Holstein Friesian cows. DNA analysis revealed 2 nonsynonymous SNPs (5174T>C, 5251C>T) on exon 4 of β-LG gene of Holstein Friesian cows, that caused a change in amino acids Isoleucine to valine and glycine to aspartic acid respectively which led to a change in protein structure and function. It was noted that genotypes TT and CT were found to have association with milk yield per day, lactose % and fat % on SNPs 5174T>C and 5099T>C, with genotype TT contributing more on milk yield per day and lactose % and genotype CT contributing more on fat %. Milk yield per day, lactose % and fat % were associated with genotypes GG and AG of SNP 4982G>A, with GG having the higher effect on milk yield per day and lactose %, while AG had great impact on fat%. CC and CG genotypes of SNP 5123C>G were associated with milk yield per 30 days, SNF, lactose % and fat %, with CC contributing more on milk yield per 30 days and SNF, whereas CG affected lactose % and fat % more. There was an association between CC and TC genotypes of SNP 5251C>T and fat %, SNF and lactose %, with CC having high impact on fat % and TC having high effect on SNF and lactose %.

The South African National Research Foundation provided financial assistance for this work (Postgraduate: Reference No. MND210614611303). The authors express their gratitude to the University for providing workspace for this research and to the Limpopo dairy farm for allowing them to gather data.
 
Author’s contributions
 
TCM, KM, TLT and VM designed the experiment, TCM carried out the data analysis. TLT, TCM, KM and VM conducted the fieldwork and compiled the manuscript. TLT updated and edited the compiled manuscript. The final manuscript was approved by all authors.

The authors declares that there is no conflict of interest.