Introduction

Domestic buffaloes (Bubalus bubalis) are of Asian origin and are known for their docility, longevity, hardiness, and adaptability in different regions of the world with distinct soil and climate characteristics. Thus, they are found inhabiting areas where cattle would be unable to be raised (flooded areas). In addition, they have high feed conversion rates, generating products of high nutritional quality through a diet with low dietary attributes (1,2). Therefore, raising dairy buffaloes is considered an agricultural activity of great financial importance, capable of providing milk with excellent nutritional properties. Initially, buffalo farming in the country was focused on beef, but from the 1980s and 1990s onwards, the scenario began to change, mainly due to the enthusiasm of producers in investigating milk production, a fact that was noticed by the growth of industrial units dedicated to the production of buffalo milk derivatives (3). This commitment was also due to the composition that buffalo milk presents compared to that of other ruminants, mainly about its composition and characteristics, especially when it comes to the levels of protein, fat, lactose, mineral residue (calcium and phosphorus), and total solids since it allows high industrial yield in the manufacture of derivatives, which can exceed 40-50% of the industrial yield when compared to bovine milk (4,5). In addition, buffalo milk is more digestible for people with lactose intolerance (6). However, even though the species displays its full potential for the activity, its production remains behind that of bovine milk. Because of this, the search for improving buffalo milk production has emerged among livestock farmers and researchers in the field, shown mainly in recent years by the growth of studies aimed at detecting variations at the level of deoxyribonucleic acid (DNA) in the buffalo genome that is capable of explaining the differences in their milk production. Among the technologies used in the search for candidate genes, chips with a high density of SNP markers and different GWAS methodologies have revealed the most genes for milk production in buffaloes (7-9). One of these is the alpha-lactalbumin gene (Lalba). It encodes whey protein (α-LA), which is essential for the synthesis of lactose in the mammary gland, which in turn is an important osmotic factor present in milk, playing a critical role in regulating the volume produced. It is also classified as an acute phase protein (APP). Since the concentration of α-LA during infections was reduced, it was then classified as APP negative (10). In a study with Polish Holstein cows, it was discovered that the SNP g.-1001T>C in the promoter region of the gene influenced milk production characteristics. Among them, high daily milk and dry matter production, in addition to high lactose production and concentration, are associated with the TT genotype (p≤0.05) (11). In Bhadawari buffaloes, the AB and BC genotypes were found to be significantly associated (p≤0.05) with total milk production and daily milk production (12). In Nilli Ravi buffaloes, 5 SNPs were detected; however, only Lalba2 at position 34310940 was associated with milk protein (p≤0.05) (5).

Although few studies have demonstrated the presence of SNPs in the Lalba gene of buffaloes, it is considered polymorphic. Because of this, researchers are increasingly using sequencing tools in different populations of the same species to prove the existence of these genes associated with zootechnical indices for milk production, which are responsible for ensuring an improvement in herd productivity. Therefore, the present research aimed to investigate and detect new possible SNPs in the Lalba gene to associate them with zootechnical data on milk production in a herd of crossbred buffaloes from the Amazon region.

Material and methods

Ethical aspects

The study was approved by the Animal Use Ethics Committee (CEUA) of the Federal Rural University of the Amazon (2820150222/2022/CEUA/UFRA).

Animals

A total of 104 crossbred Murrah and Mediterranean buffaloes aged between 1 and 5 years weighing between 450 and 500 kg were used. All were healthy and free of mastitis and belonged to a property located in the municipality of Bujaru, Belém, Pará, Brazil (1°37’40’’S, 48°13’21’’W). Five ml of blood were collected from these animals in tubes containing anticoagulant (EDTA).

For the evaluation of gene expression, 42 buffaloes with average milk production ranging from 2 to 10 L/day were selected. From these, 10 mL of milk were collected in sterile 15 mL Falcon tubes, free of DNase, RNase, and non-pyrogenic, with 2 mL of RNAlatter added to stabilize and protect the RNA in the samples. Both samples were stored in a thermal box at a maximum temperature of 4°C and taken to the Serology and Molecular Biology Laboratory of the Federal Rural University of the Amazon – UFRA for further analysis.

DNA Extraction and Polymerase Chain Reaction (PCR)

It was performed using the phenol-chloroform method (13) and stored at 20 °C. The investigation of genetic variation based on SNPs was performed through PCR of the Lalba gene (ID: 102410146) in its promoter region and 5’ UTR. The following primer pair F 5’ ACC TTT CCC AAA AAG GTT GG 3’ and R 5’ ACC ACA CTG CTC ACC TCC TT 3’ used in the research was designed through Primer3plus based on the reference sequence deposited in Genbank, amplifies a fragment of 833 bp. The reactions were performed for a final volume of 25 μL, being 2 μL of DNA, 12.5 μL Taq DNA Polymerase 2x Master Mix (Ampliqon, Odense, Denmark), 1.0 μL of each primer - forward and reverse - (10 pmol/μL) and 8.5 μL of ultra-pure water. All reactions were performed in a Thermal Cycler 2720 thermocycler (Applied Biosystems, Foster City, CA, USA), with the following amplification protocol: initial denaturation at 95 ºC for 5 min, followed by 35 cycles with denaturation at 95 ºC for 30 s, annealing at 55 ºC for 30 s, extension at 72 ºC for 50 s and final extension at 72 ºC for 5 min. The PCR products were analyzed using 0.8% agarose gel.

DNA Purification and Sequencing

The PCR products from 83 samples were purified using the PCR Product Purification Kit on a Column (Ludwig Biotec LTDA, Alvorada, RS, BR) and visualized on 0.8% agarose gels. These purified products were sequenced using the Sanger method using the BigDye Terminator v3.1 Cycle Sequencing Kit (Invitrogen California/USA) in a final volume of 10 μL in an ABI 3500 Genetic Analyzer automatic sequencer (Applied Biosystems). The sequences obtained were edited using the Finch TV software version 1.4.0 (Geospiza Research Team, USA), compared with the reference sequence for buffaloes entered in GenBank in the BLAST system, and finally aligned in BioEdit 7.2 using the ClustalW tool. Transcription factor domains in the promoter region were identified using Nsite v. 5 software (14).

RNA Extraction

For RNA isolation, milk samples were processed to separate somatic cells. Then, 1000 μL of milk from 15 mL Falcon tubes added with RNAlatter were centrifuged at 12,000 rpm at 4°C for 8 minutes. The supernatant was carefully removed with a spatula. The remaining material was transferred to sterile 2.0 mL Eppendorf tubes, where 700 μL of Trizol® (Invitrogen, USA) was added. The extraction protocol was followed according to the manufacturer's recommendations. The samples were quantified in a BioDrop μLITE spectrophotometer (BioDrop) using the A260nm/A280nm ratio.

Gene expression (RT-qPCR)

The gene expression profile was investigated from the 42 buffaloes in exon 2 of the Lalba gene. For this, the primer pairs described in table 1 were designed using Primer3Plus based on the reference sequences deposited in GenBank. The RT-qPCR reactions were performed in triplicate for Lalba and the endogenous Glyceraldehyde-3-phosphate dehydrogenase (GAPDH) gene. The amplification reaction protocol was performed for a final volume of 10 μL according to the recommendations of the manufacturer of the RealQ Plus and RealQ Fast kit (Ampliqon, Odense, Denmark). Relative expression values ​​were determined by the equation 2-∆Ct, where ∆Ct is the ratio between the Ct of the target gene and the Ct of the endogenous gene (15).

Table 1: Primer pairs designed in Primer3plus according to the reference genes deposited in GenBank (ID) for gene expression analysis in dairy buffaloes

Statistical analyses

Genetic diversity data such as allele and genotypic frequencies, Hardy Weinberg Equilibrium (HWE), and F_IS were obtained using PopGene (16). Expression patterns were analyzed using SAS OnDemand, assuming a normal distribution. Mean milk production data were analyzed using simple descriptive statistics. The Tukey test was used to assess the association with polymorphism data and the gene expression profile of haplotypes, with a significance level of 0.05.

Results

The alignment performed on all the sequences obtained revealed the presence of 10 SNPs, seven of which were transitional and two of which were transversional. In addition, the presence of an insertion block with the nucleotide sequence TAAA between positions -552 and -555 was also detected in 92.8% of the samples (Table 2).

Table 2: Summary of single nucleotide alterations (SNPs) in the LALBA gene in buffaloes from the Bujaru-Pa region compared to Bubalus bubalis (GenBank ID: 102410146)

UTR = untranslated region; bp: base pairs

The analysis of population genetic diversity detected that only eight SNPs presented the frequency of the wild-type allele greater than 0.5, equivalent to 90% of the total. These alleles of the SNPs in positions -339 (97.6%), -60 (92.5%), and -720 (91.5%) were those that exhibited high frequencies in the population that exceeded 90%. Regarding the frequency of the mutant allele, SNPs -42, -444, -501, and -591 were the exceptions, as they were the only ones that exhibited the frequency of the mutant allele (G) greater than 0.5, being 68.1%, 95.1%, 95.1%, and 97.5%, respectively. Furthermore, the insertion mutation was detected in 77 buffaloes evaluated between positions -555 to -552, which contributed 4 new nucleotides (TAAA) to the population, also at a high frequency of 96.3%. Six of the 83 buffaloes are heterozygous for this SNP, that is, characterized as InDel, and no buffalo in the present study presented the deletion, as verified in the reference sequence for the species deposited in the GenBank database. All values ​​for the effective number of alleles in the population for all SNPs evaluated were above 1.

The SNPs -42, -292, -445, and -720 presented lower observed heterozygosity rates (Ho) than expected (He), according to the HWE proportions, which is evidenced by the positive F_IS values, suggesting an excess of homozygotes, with the inbreeding process occurring for these. The only SNPs that exhibited the presence of homozygotes for the mutant allele were -42, -292, -445, and -720. On the other hand, for loci -60, -339, -444, -501, -591 and the insertion block together, the Ho rates were higher than those of He, data confirmed by negative F_IS values, thus suggesting an excess of heterozygotes and possible selection of them in the population. Locus -591 presented the lowest value of 0.037. This analysis of the SNPs found suggests the presence of low diversity. However, loci -42, -292, and -445 presented the highest values ​​of 0.626, 0.601, and 0.601, respectively, which may suggest the presence of greater diversity for these loci within the population. Regarding the proportions of HWE, of all the SNPs found in the present study, nine are under equilibrium conditions (P>0.05), which are in the following positions -42 (A>G), -60 (T>C), -292 (G>A), -339 (C>T), -444 (T>G), -445 (C>T), -501 (T>G), -555 (TAAA) and -591 (C>T). The SNP in position -720 (C>T) is deviating from the HWE (P<0.05). All the information mentioned above is described in table 3.

From the analysis of the binding to transcription factors performed by the software above, it was possible to see that of the ten SNPs found, only one, at position -42 (TCATAAATA), has a strong binding site for the transcription factor hepatocyte nuclear factor 3 (HNF3) in its promoter region. However, the same analysis also found five other transcription factors, namely, Myocyte enhancer factor 2 (MEF2), Specificity protein 1 and specificity protein 3 (SP1; SP3), Sterol regulatory element-binding protein (SREBP), and Zinc finger proteins and Myc-associated zinc finger protein (MAZ/Sp1). However, these are not associated with any of the nine SNPs but are binding sites for the regulatory region of the Lalba gene and are, therefore, important to mention.

Table 3: SNP positions and respective descriptions of genetic diversity in the promoter region and 5' UTR of the Lalba gene in dairy buffaloes

Ho: Observed heterozygosity, He: Expected heterozygosity, Ne*: Effective number of alleles, I*: Shannon diversity index, F_IS: Inbreeding coefficient, HWE: Hardy Weinberg equilibrium (P>0.05)

The association of the SNPs found with the average milk production (L/day) of the buffaloes evaluated is described in table 4. The genotypic comparisons of all the SNPs found in the study with the average milk production did not show any statistical difference (P>0.05). However, even though it was not possible to observe a significant effect in the association, it was possible to notice that some SNPs stood out because they were responsible for increasing milk production in certain buffaloes.

Therefore, for SNPs at positions -60, -339, -501, -720 and for Indel (-555), it was possible to observe that all heterozygotes of the same exhibited greater milk production compared to their wild homozygotes. Highlighting the heterozygote (CT) of SNP -720, which presented the highest milk production observed among all other SNPs. In addition, when compared to the production exhibited by its wild homozygote (CC), there was an increase in production of more than 2 liters of milk for the heterozygote, being 7.44 L/day; the homozygote (TT) also followed this increase, as it exhibited 7.01 L/day. Based on these findings, we can infer that the presence of mutant alleles in these SNPs may be contributing positively to the increase in milk production, having the so-called additive effect.

This deduction can be better observed and explained in SNPs -444, -501, and -591, for which it was not possible to find any buffalo belonging to the respective wild homozygous genotypes. However, the mutant homozygotes all demonstrated an increase of more than 0.5 liters in milk production compared to their heterozygotes. SNPs -42, -292, and -445 were the only ones that presented the three genotypes in the population, with the buffaloes of the wild homozygous genotypes exhibiting higher milk production.

Table 4: Association of novel polymorphisms detected in the regulatory region of the Lalba gene with buffaloes' average milk production (L/day)

Regarding the analysis of the gene expression profile performed in the present study, all haplotypes evaluated exhibited expression of the Lalba gene, as well as the selected endogenous gene, GAPDH, in the somatic cells of milk. Haplotypes three (ATGCGCGINTC), four (GTCGACGCTGINTC), and two (AGTGCGCGINTC) presented the lowest (0.632±0.263) and highest (18.986±16.264; 15.386±12.956) expression levels, respectively. The last two mentioned were statistically significant (P<0.05) when compared to haplotypes one, three, five, and seven. The only one that was similar to the significant ones was haplotype six, which exhibited an expression profile of 9.706±5.199, a value close to those with the highest expression (Figure 1).

The gene expression profiles obtained from the seven haplotypes were compared with the average production data of the buffaloes evaluated in the study. However, as with the ten new SNPs found, there was no statistically significant difference in the average milk production (P>0.05). However, based on the analysis performed, it was possible to perceive that haplotype four, which presented the highest expression, could be considered a candidate for greater synthesis of mRNA of the Lalba gene in milk and, consequently, greater production of the α-lactalbumin protein in buffaloes.

Figure 1: Relative gene expression from mRNA of the Lalba gene obtained from somatic cells of milk of the seven different haplotypes. According to the equation 2-ΔCt. Different letters above the bars indicate significant pairwise comparisons according to the Tukey test (P<0.05).

Discussion

The characteristics related to milk production in ruminants are controlled by quantitative loci traits (QTLs). These are spread throughout the genome and are regions in which genotypes are associated with phenotypes demonstrating differences in milk production and composition, for example. Therefore, the wide variety of cows selected for greater and better production and composition, both in quantity and quality, depends on identifying several candidate genes and their polymorphic variants (17,18). The search for so-called promising genes is commonly used due to the ease of understanding certain metabolic pathways and physiological processes and because they are close to QLTs (5).

The Lalba gene has been reported as polymorphic in many production animals, such as goats, sheep, cattle, and buffaloes. Therefore, there is evidence that nucleotide substitution mutations in regulatory and non-coding regions of genes influence the definition of complex characteristics (5,11,19,20).

In the present study, it was possible to perceive the variety of SNPs that can be found in buffaloes since ten new SNPs were detected in the regulatory region of the Lalba gene, with the one at position -42 (A>G) presenting a binding site for the transcription factor hepatocyte nuclear factor 3 (HNF-3). HNF-3 is a factor that plays a crucial role in the regulation of metabolism and tissue differentiation, mainly in the pancreas and liver (21). In another study with dairy cattle, it was observed that the SNP at chromosomal position g.146704699 G>A of the AGPAT3 gene, especially the A allele, created a transcription binding site (TFBS) for FOXA1 (Forkhead box protein A1, also known as hepatocyte nuclear factor 3-alpha (HNF-3-alpha). In addition, it showed a significant effect (P<0.05) on three short-chain fatty acids (C6:0: caproic fatty acid, C8:0: caprylic fatty acid, and C10:0: capric fatty acid) present in the milk of Holstein cows (22).

Regarding the genetic diversity parameters evaluated, the values ​​found for allele and genotypic frequencies of wild-type alleles in most SNPs were high, which supports the degree of consanguinity observed in some SNPs through positive F_IS values. However, even with these SNPs suggesting an inbreeding process, 90% of those evaluated in the study are under HWE conditions due to the proximity in He and Ho values. Furthermore, this fact can also be explained by the assumption that inbreeding does not necessarily eliminate all alleles from a population but rather reorganizes them when present. Significant differences in He and Ho were observed only for the SNP (-720 C>T), which presented a deviation in HWE (P>0.05), so that the population may be subdivided by the presence of significant inbreeding or by the displacement of genes even from another population (23,24).

Many studies aimed at verifying genetic variability in water buffaloes are being carried out in many regions of the world; however, most are focused on microsatellite markers. In a study carried out with seventeen populations of Turkish water buffaloes, it was found that all twenty microsatellite loci tested were highly informative and polymorphic. However, the average F_IS value obtained in the study was 0.091, with all loci in HWE deviation (25). In the present work, the average F_IS obtained was 0.098, similar to what occurred in the previous study, but with divergence regarding HWE where of the ten SNPs found, nine are under HWE and only one in deviation. However, regarding the F_IS values ​​isolated from each of the SNPs found, it can be suggested that the population studied presents only four in the process of inbreeding (positive values), characterizing a reduction in diversity based on many causes, such as the Wahlund effect, which consists of an increase in homozygotes and a decrease in heterozygotes in samples of individuals belonging to subdivided populations (26).

Regarding genetic diversity studies carried out with buffaloes in the Amazon region, the F_IS values ​​found in the present work were higher than those found by other researchers in the literature (27-29), 0.0294 to 0.0329, 0.0741, and 0.0730 respectively. To investigate the genetic variability of the SNP (LEP-1620) in the leptin gene, a study was carried out with two buffalo breeds (Mediterranean and Murrah) and their crossbreds in two distinct breeding systems characteristics of the Amazon region (dry land and floodplain). In these populations, distinct F_IS values ​​were found, being negative (-0.125, -0.171, and -0.142) for Mediterranean and Murrah from dry land Mediterranean from floodplain and positive (0.057) for crossbreds raised in a floodplain system with all racial groups in HWE (P>0.05) (30). And the leptin gene is strongly related to milk production and characteristics such as fat and protein in buffaloes (31).

Research based on SNP detection in the regulatory region of the Lalba gene has been studied since the 1990s, where many SNPs have been reported since then in several ruminant species such as cattle, goats, sheep, and buffaloes, most of which are significantly associated with increased milk production, in addition to being associated with the main components of milk such as proteins and fats. And because they are located in the regulatory region, they can also influence the binding to many transcription factors and the levels of gene expression profiles (11,32-36).

Association studies with milk production traits for the Lalba gene are still scarce; however, one of the first studies aimed at elucidating the effect of an SNP in exon 1 of the gene in two distinct buffalo breeds, Bhadawari and Murrah, was able to demonstrate that the AB and BC genotypes present in Bhadawari individuals were significantly associated (P<0.05) with total milk production and daily milk production. Still, the AB, BB, BC, CC, and CD genotypes present in Murrah did not show significant effects on these traits (12). Another non-significant SNP, the α-LA2516 SNP, located in exon 4 of the gene, was found to have no relationship with protein percentage, fat percentage, milk production, and number of somatic cells in milk (P>0.05) in Chinese Holstein cows (37).

Other authors have also found SNPs in the Lalba regulatory region, both in its 5' flanking region and in intron regions, which also had no significant effect on performance and milk production traits in cattle (38,39). This corroborates the findings of the present study, in which no significant association was observed regarding average milk production (L/day). The lack of significance in the associations with the animals' dairy performance may be due to the environmental, epigenetic effect, or intrinsic differences that the studied breeds have, which, in the case of buffaloes, are still little investigated (5).

However, more recent studies have demonstrated that there are significant effects between associations of certain types of mutations such as SNPs, copy number variants (CNVs), and even microsatellite studies in both river and swamp buffaloes (5,25,40). Regarding environmental and epigenetic effects, variations based on CNVs contribute to better environmental adaptation in buffaloes across the continent where they are raised. Thus, it was reported in the literature that a CNV was observed in the promoter region 1647 bp upstream of Lalba, and this overlaps with two different transposons (RTE-BovB and Bov-Ta3). But the interesting thing about this finding is that this CNV was common (CN = 2) in 95.08% of river buffaloes but was presented as a deletion (CN = 0) in 82.50% of swamp buffaloes, which may have contributed to the higher milk production of river buffaloes compared to swamp buffaloes.

SNP mutations can be found throughout the genome of a wide variety of species. However, when they are located along the regulatory region, they can interfere with the binding of existing transcription factors to the target DNA sequence analyzed, thus modulating the gene expression process through transcription inhibition or activation (41). Concerning animal production, this phenomenon can result in divergent reproductive and productive efficiencies in milk production.

The isolation of epithelial cells derived from somatic milk cells (MECs) is being used as a suggestion for non-invasive collection methods to obtain expression analysis of genes involved in milk production as well as its components. High amounts of mRNA from protein and fat genes were observed in MECs obtained from the Murrah breed compared to the Sahiwal cattle breed. The expressions for casein and serum protein (Lalba) genes were significantly (P<0.05) higher in MSCs purified during the early stage of lactation compared to those purified from other stages (peak, middle, and late) (42).

Another alternative to obtain mRNA of milk production genes is the somatic cells present in milk (MSCs), which is done by separating these cells by centrifugation. This was the method used in the present study to obtain mRNA from MSCs in the milk of buffaloes selected for evaluation of the expression of the Lalba gene. In the study carried out in buffaloes at the end of lactation, the expression of Lalba was low in MSCs; however, it was higher and significant (P<0.05) in buffalo milk fat globules (MFG), indicating that MFG can be used as another alternative for collecting these transcripts since the total RNA present in MFG derives from mammary epithelial cells (43).

From transcriptome analyses of four dairy buffaloes, it was possible to demonstrate that genes involved in production were highly expressed during all phases of lactation (initial, middle, and final). One of them was the Lalba gene, which exhibited high levels of expression in early lactation, with a subtle decline in the middle stage, but maintained its increased expression at the end (44). The increased expression of this gene is relevant to the synthesis of milk proteins since α-lactalbumin generally makes up about 20% of whey proteins (45). In addition, α-lactalbumin encoded by the Lalba gene is responsible for regulating lactose synthesis in many mammals, so lactose is one of the most valuable carbohydrates present in milk and plays an important role in the regulation of osmolarity. Therefore, the Lalba gene is considered relevant to the lactose production pathway in milk (46,47).

The data shown in the literature corroborate the gene expression values ​​of the Lalba transcripts obtained in the present study, which were derived from milk MSCs, and all the haplotypes evaluated showed expression of the same, also showing that MSCs are excellent sources that provide transcripts to be evaluated in gene expression analyses. In which of the seven haplotypes evaluated, sets four and two were highly expressed and significant (P<0.05) in the MSCs of the crossbred buffaloes evaluated with levels much higher than those demonstrated by the authors cited above.

Conclusion

This study was able to demonstrate that the Lalba gene is polymorphic in crossbred dairy buffaloes raised in the Amazon region. Although no significant association with milk production was found, these haplotypes may have contributed to the significant expression of serum protein mRNA in them. Lalba gene transcripts in MSCs are expressed in buffalo milk. Of all the SNPs found, InDel stands out as an important piece of information for the genetic characterization of the population. This is present in 7% of the population and, in addition, approximately 93% of the animals are characterized by the presence of insertion of four nucleotide bases (TAAA) in their regulatory region, which can contribute with great importance to the evolutionary processes of living beings, as it can guarantee genetic variability in addition to the emergence of new adaptive characteristics. However, more associative studies should be carried out on buffalo breeds, especially those raised in the Amazon region, to further promote and improve the dairy production chain.

Conflict of interest

This manuscript's authors declare no conflict of interest regarding the writing process or data analysis.

Acknowledgments

Gratitude to the team of the Molecular Biology laboratory of the Federal Rural University of Amazon and the financial support of this research by the Coordination for the Improvement of Higher Education Personnel Foundation (CAPES).