CSN1S1, CSN3 and LPL, Three Validated Gene’s Polymorphisms Useful for a More Sustainable Dairy Production in the Mediterranean River Buffalo

1. Introduction

The domestic water buffalo (Bubalus bubalis) is a tropical animal characterised by a marked ability to adapt to the environment and a high efficiency of feed use in conditions of forage shortage. The species originated in Southern East Asia where nowadays 97% of world buffaloes in world still are reared [1], and spread west arriving in Syria, Egypt and then west Europe [2]. Therefore, these animals are of major economic and cultural importance for many populations worldwide, supplying milk, meat and draught power. Two buffalo sub-types exist, the swamp type (2n=48) exclusively present in the native Asian continent and the river type (2n=50) globally more spread also in the other continents. These buffalo sub-types differentiate for karyological, morphological and behavioural characteristics [3,4,5].

Italy is the European country with the greatest number of buffaloes raised. In recent years, the Italian buffalo population has increased from about 12,500 heads in the 1950s to over 400,000 in 2019 [1] that represent about 85% of the entire European population. Such expansion took place thanks to the exploitation of the buffalo milk by the national and international increase in the “Mozzarella di Bufala Campana PDO” demand. Recent data show a significant growth of the whole supply chain with a turnover estimated at 500 million euros, more than 20,000 operators and a 5% annual export increment (www.ismea.it). Despite that, compared to other ruminants, the domestic buffalo have received less attention and economic investments; therefore, the species possess a great improvement potential.

The achievement of high production levels and good efficiency implies the optimization of a number of factors and processes including genetic improvement. In this respect, although new knowledge has been acquired [6,7], a high contiguity assembly of the reference genome has been published [8], and the first SNP array designed specifically for buffaloes has become available [9], the use of genomic data is still very limited. Therefore, nowadays, the estimation of genomic breeding values and the application of genomic selection have huge delays in domestic buffalo, as recently reported also by Cesarani, et al. [10]. In addition, the genome-wide (GWAS) approach in buffalo by the medium density 90K array SNP often identifies candidate variants in intergenic regions nearby many potential genes of interest [11,12,13], but no further confirmation studies are then carried out. For this reason, the candidate gene approach is still today a valid method for the identification of genetic associations with milk production traits. At the same moment, it is a useful information for the associations of breeders that, in the last decade, promoted the selection of buffalo sires with favorable genotypes for milk traits (https://www.risbufala.it/?page_id=58841).

Milk yield [14,15,16,17,18], total protein and caseins [19,20,21,22,23], fat content, fat percentage and fatty acid composition [18,20,24,25,26,27,28,29,30], milking time [14] etc., are among the most studied traits and those of great attention for the breeders’ associations because directly connected to cheese yield and economic profit. Genetic variability and association with dairy traits have been found for many genes of economic interest (CSN1S1, CSN1S2, CSN3, SCD, LPL, OXT, OXTR, etc.). However, many association studies are often investigated into single buffalo farms, with a limited number of samples, or carried out using single genes’ variants. For instance, the association between the protein percentage and AJ005430: c.578C>T on CSN1S1 (αs-1 casein) [21] or the milk yield and the FM876222: g.133A>C on SCD (Stearoyl-CoA Desaturase) [15]. Therefore, aim of this study was to extend the genotyping of the most four promising SNPs in 4 genes of interest for selection goals (CSN1S1, CSN3, SCD, LPL) in a larger population and validate genetic relationship with milk traits for breeding purpose.

2. Materials and Methods

2.1. Sampling and DNA Isolation

Individual blood collection was performed in compliance with Italian national laws and regulations by official veterinarians of ASL (Local Sanitary Unit of the Ministry of Health) during the routine farm prophylaxis.

Sample collection was carried out on a total of 800 Italian Mediterranean river buffaloes belonging to 8 dairy farms mainly located in Campania region (Southern Italy). Individual blood samples were collected during the routine farm prophylaxis by official veterinarians of ASL (Local Sanitary Unit) of the Ministry of Health.

Genomic DNA was isolated using the procedure described by Goossens and Kan [31]. Concentrations and OD_260/280 ratios were measured with the Nanodrop ND-2000 Spectrophotometer (Thermo Fisher Scientific Inc., Waltham, MA, USA).

2.2. Genotyping

Genotyping was accomplished using PCR based methods described by Pauciullo et al. [22] for the AJ005430: c.578C>T at the CSN1S1 (αs1-CN) and HQ677596: c.536C>T at the CSN3 (k-CN); Gu et al. [28] for the FM876222: g.133A>C at SCD (Stearoyl-CoA desaturase) and Gu et al. [27] for the AWWX01438720.1: g14229A>G at LPL (Lipoprotein lipase). PCR amplification was carried out using BioRad T100 thermocycler (BioRad). The digestion products were analysed directly by electrophoresis in 2.5% agarose gel in 0.5× TBE buffer and stained with SYBR green nucleic acid stain (Lonza Rockland Inc.) (Figure S1).

2.3. Phenotypes Collection and Dataset Editing

The phenotypes for milk yield and composition from official recording program of the Italian Association of Breeders (AIA) were used for this study under a cooperation agreement. Daily milk yield (dMY, kg), protein yield (dPY, kg) and fat yields (dFY, kg), fat and protein contents (dFP, % and dPP, %), somatic cell count (SCC, 10³cell/ml) and urea (mg/dl) were individually recorded in each month for the whole lactation from 2010 to 2021 (n=16,457 records). Only animals with both complete genotypes for the 4 SNPs (MAF>0.05, n=762) and lactation records were retained as valid records in successive analysis. A total of 15,742 individual milk test day records (2,496 lactations) were available on 680 buffalo cows with 3.6±1.7 parity (from 1 to 13) and 6.1±1.2 test day records per lactation on average (Table 1).

A further data editing has been performed prior to study the association between SNP genotypes and milk phenotypes: i) to remove outlier unsounded data (greater than ±3.5 standard deviations or null values per DIM>10 d), ii) to transform SCC in SCS according to Ali and Shook [32] and iii) keeping those buffalo cows whose lactations had at least 5 records. The final dataset included 645 buffaloes.

2.4. Statistical Analyses

Descriptive statistics were performed both on SNP and phenotypic data. Minor allele frequency and Hardy-Weinberg equilibrium test were computed for all 4 genes. The pairwise Pearson correlations among milk traits were also calculated.

The effect of SNP genotypes on milk traits were assessed through mixed model analysis implementing 2 different genetic model using proc mixed of SAS ^® (2016 Cary, NC, USA): allelic and genotypic models. In both genetic models, the fixed effects of contemporary group were included as well as other systematic source of variation, as specified thereafter. The type 3 sum of square of proc mixed were computed and SNP effects were considered significant for p-values <0.05.

2.5. Allelic Model

In the first genetic model for each of the 4 investigated polymorphisms, the phenotypic values for milk traits were regressed onto the number of B allele (0, 1 or 2) for A/A, A/B and B/B genotypes respectively (Table 2) according to an additive model. Moreover, the effect of dominance was assessed with a different parametrization for dominant (1) and recessive genotypes (0). The general model used both for additive and dominance parametrization was:

$$y_{ijklmo}=\mathsf{\mu }+\mathsf{\beta }\times \mathrm{S}\mathrm{N}\mathrm{P}+{\mathrm{Y}\mathrm{E}\mathrm{A}\mathrm{R}}_i+{\mathrm{D}\mathrm{I}\mathrm{M}}_j+{\mathrm{N}\mathrm{L}}_k+ {\mathrm{S}\mathrm{E}\mathrm{A}}_l+{htd}_m+{bcow}_n+e_{ijklmn}$$

(1)

where y is the test-day phenotypic values for each analysed milk trait, SNP is the covariate of allelic count and $\mathsf{\beta }$ the average substitution effect for the additive model (AM), or dominance effect for the dominance model (DM). Moreover, the fixed effect of year of birth (11 levels), days in milk (DIM: 15 classes of 20 d each), parity (NL, six classes, from 1 to 6+) and season of birth (SEA, 2 classes: autumn-winter and spring-summer) were fitted. Random effects for combination of herd-test day (htd, 669 levels), buffalo cows (bcow, 646 levels) and residual were also included. Random effects were assumed independently and identically distributed.

2.6. Genotypic Model

The genotypic was alike the model (1) with the main difference that the genotypes at the four loci were treated as cross classified fixed effect instead of covariate (3 genotypic class for A/A, A/B, B/B) according to:

$$y_{ijklmo}=\mathsf{\mu }+{\mathrm{S}\mathrm{N}\mathrm{P}}_i+{\mathrm{Y}\mathrm{E}\mathrm{A}\mathrm{R}}_j+{\mathrm{D}\mathrm{I}\mathrm{M}}_k+{\mathrm{D}\mathrm{I}\mathrm{M}\left(\mathrm{S}\mathrm{N}\mathrm{P}\right)}_{k\left(i\right)}+{\mathrm{N}\mathrm{L}}_l+{\mathrm{S}\mathrm{N}\mathrm{P}\left(\mathrm{N}\mathrm{L}\right)}_{i\left(l\right)}+ {\mathrm{S}\mathrm{E}\mathrm{A}}_m+ {htd}_n+{bcow}_o+e_{ijklmno}$$

(2)

where DIM(SNP) and SNP(NL) are the nested effect of days in milk within SNP genotype and the genotypes nested within parity effect, and the other are the same as above. In this model and type 3 sum of square F-test for fixed effect were performed and the marginal means of different genotype are separated at p-values <0.05 in post-hoc comparison adjusting the p-values according to Tukey HSD (adjust=Tukey of proc mixed).

Finally, to estimate the proportion of variance explained by genotypes a simplified model was used:

$$y_{ijklmo}=\mathsf{\mu }+{\mathrm{Y}\mathrm{e}\mathrm{a}\mathrm{r}}_i+{\mathrm{D}\mathrm{I}\mathrm{M}}_j+{\mathrm{N}\mathrm{L}}_k+ {\mathrm{S}\mathrm{e}\mathrm{a}\mathrm{s}\mathrm{o}\mathrm{n}}_l+ {SNP}_m+{htd}_n+{bcow}_o+e_{ijklmno}$$

(3)

where SNP genotypes, htd and cow are treated as random effects and the proportion of variance explained by the SNP genotype ($r_{SNP}^2$), herd-test day ($r_{htd}^2$) and buffalo cows ($r_{bcow}^2$) were computed respectively as ratio of the components to the total variance for each polymorphism

$$\begin{array}{l}{\widehat{\sigma }}^2={\widehat{\sigma }}_{SNP}^2+{\widehat{\sigma }}_{htd}^2+{\widehat{\sigma }}_{bcow}^2+{\widehat{\sigma }}_e^2:\\ r_{SNP}^2={\widehat{\sigma }}_{SNP}^2∕{\widehat{\sigma }}^2,\\ r_{htd}^2={\widehat{\sigma }}_{htd}^2∕{\widehat{\sigma }}^2 \mathrm{and}\\ r_{Bcow}^2={\widehat{\sigma }}_{bcow}^2∕{\widehat{\sigma }}^2.\end{array}$$

(3)

3. Results and Discussion

In the present study, four SNPs (AJ005430:c.578C>T, HQ677596:c.536C>T, FM876222:g.133A>C and AWWX01438720.1:g14229A>G) each in a gene of interest for selection goals (CSN1S1, CSN3, SCD and LPL respectively) were genotyped in a Mediterranean river buffalo population of 800 animals belonging to 8 farms (Figure S1). The specific choice of these SNPs was driven by the need to confirm their impact on milk traits carried out in previous studies on relatively small buffalo populations [15,21,23,27]. In addition, two of these SNPs (AJ005430:c.578C>T in the CSN1S1 and HQ677596:c.536C>T in the CSN3) were recently included in the genotyping program of buffalo sire selection by one of the two Mediterranean buffalo associations of breeders (Research Innovation and Selection for the buffalo).

The four investigated SNPs largely segregate in the buffaloes’ population under study (MAF >0.21, Table 2) with a range of variability of 0.16-0.55 across genes or herds. With few exceptions, the four polymorphisms were in HW equilibrium within and across herd (Figure S2, Table S1). Overall, the deviation from the HW equilibrium was partially expected for the SCD (χ2=6.19) that was previously investigated in two different populations with similar findings (χ2=6.92, [15]; χ2=7.96, [28]). The SCD FM876222: g.133A>C was associated with milk yield and the allele substitution effect was assessed in about -1kg/d with 12% of the total phenotypic variance explained by the polymorphism [15]. Such an effect is larger than that evidenced for the DGAT1 on milk yield in dairy cattle [33]. Despite that, so far, no marker assisted selection was voluntarily applied in favour of the allele A to increase the buffalo milk production. Therefore, the HW deviation for the SCD, with the frequency of the allele A almost reaching 80%, can be considered as the result of farmers’ directional selection for more productive animals.

Conversely, the deviation of the CSN1S1 from the HW principle (χ2=5.06) was unexpected considering the previous studies [21,22]. However, starting from 2021, the Italian buffalo population is under selective pressure for the SNP AJ005430:c.578C>T (https://www.risbufala.it/?page_id=58841). Therefore, potentially, the HW deviation can be considered as result of an artificial selection sweep.

For six milk traits the number buffaloes with valid records were 646 with an average DIM of 153±93 d, whereas 29 animal were discarded for SCS due to missing phenotypes. The number of test days and lactations records slightly differs for the milk traits (from 20 to 22 records per animal on average). The average daily milk yield and composition and their pairwise phenotypic correlations (Table 3) are in accordance with previous reports [10,14,15,19,34,35] and with the official milk yield mean (8.70 ± 2.58 kg/d) reported for standard lactations (until 270 DIM) in 2022 [36]. Milk urea, important for its role in nitrogen metabolism, shows a weak correlation (<0.10) with all traits. Indeed, milk urea correlated positively with protein yield and negatively with fat contents (Table 3).

This result is among the first indication of correlation between milk urea and other milk parameters in buffaloes, since little studies are available in this species. Instead, more information is available in dairy cows, as well as more conflicting data are reported. In general, a low negative genetic correlation was found between milk urea and milk yield [37,38], but in New Zealand dairy cattle the correlation between these two traits was reported as moderately positive [39,40]. Differences between diet formulations are considered as important elements that may cause genetic × environmental interaction that could explain such differences [37]. This could be also the case of the buffalo, whose genetic background is different from the dairy cattle, as well as the energy requirement and diet.

With few exceptions (dFP and SCS in respect to birth season) all the fixed effect were highly significant (Table S2). Additive and dominance effect were reported in Table 4. In the Allelic models, LPL show a significant negative substitution effect on dMY when increasing the number of G alleles (p<0.05) and positive effect on milk contents of fat and protein (dFP and dPP p<0.05).

Considering that the lipoprotein lipase (LPL) facilitates the hydrolisis of triglycerides transported via chylomicrons and very low-density lipoproteins, serving as the pivotal stage in the transportation of free fatty acids to mammary gland and adipose tissues, through its regulation of fatty acid delivery to the mammary gland, LPL could influence the fat content of milk.

Our result is also consistent with the recent findings in the Italian buffalo population. In fact, the allele G in homozygosis showed a significant over expression (~2.5 fold higher) compared with other genotypes and it was associated with milk PUFA content [27]. Conversely, the allele A in homozygosis showed higher values for the milk yield, although the estimated difference with the other two genotypes only approached the level of significance (P=0.07) [27]. Associations of LPL with milk fat traits and dMY were also found in other species [41,42,43,44]. So far, no associations between LPL and milk proteins were reported for buffaloes, but recently in Czech dairy goats a significant association was found for this trait for the SNP LPL g.185G>T [42].

The investigated polymorphism at CSN1S1 exhibited positive additive effects on dMY, dFY, dPY and SCS at increasing dose of T alleles (Table 4), whereas no significant effect of CSN3 polymorphism was exerted on proteins (dPY and dPP) and other milk traits (dMY, dFY and urea), with the exception of higher SCS observed at increasing number of T allele (Table 4). Overall, this result confirms and reinforces the importance of the αs1-CN encoding gene in the determination of buffalo milk characteristics with some important differences compared to the former study of Cosenza et al. [21]. The first is the higher number of associated dairy traits found in the present study with the same SNP, although the protein percentage showed only a tendency in the genotypic model (p<0.09), whereas associated (p<0.04) by Cosenza et al. [21]. However we can consider the present dataset more robust (2500 lactations, 8 farms, nearly 650 buffaloes) compared to the former study that was numerically much lower (500 lactations, 1 farm, 175 buffaloes). This difference had also other consequences. The most important is the allele substitution effect (cytosine into thymine) that changed from -0.014 observed by Cosenza et al. [21] to 0.011 of the present study. Differences of substitution effects across populations are possible and they are function of several elements like the extent of variances (additive, dominance and additive by additive), the genetic distance of the populations and their heterozygosity [45]. The contribution of the AJ005430:c.578C>T to the total phenotypic variance found by Cosenza et al. [21] was quite low (r²_αs1 =0.003) compared to the present study (r²_αs1 = 0.100). If we further consider that Cosenza et al. [21] also found a large dominance effect (–0.028 ± 0.019), then altogether these data may explain, at least partially, the different results between the two studies.

The approached association (P<0.06) of the CSN3 (κ-CN) in the genotypic model represent a further confirmation of the importance of this locus for milk traits. The HQ677596:c.536C>T, alleles X1 (p.Ile¹³⁵) and X2 (Thr¹³⁵) are known to play a fundamental role in the buffalo milk processing, especially in combination with the variants AJ005430:c.578C>T, alleles B (p.Ser¹⁷⁸) and A (Leu¹⁷⁸) at the CSN1S1 [19,23]. In this respect, the combined genotypes AA-X1X2 showed better curd performances as shorter rennet coagulation time, shorter curd-firming time and larger curd firmness [19]. Instead the combination of the alleles CSN1S1*B and CSN3*X1 resulted in a greater curd yield [23]. Surprising was the association evidenced between both casein genes (CSN1S1 and CSN3) and SCS. The allelic and genotypic models converged in defining the polymorphism at CSN1S1 gene for both additive (p<0.05) and dominance (p<0.05) effects on SCS. The average values for C/T and T/T buffaloes did not differ each other in the log-transformed somatic cell count at p<0.05 (3.28 and 3.25) whereas the average for C/C genotypes was significantly lower than the formers, thus configuring a degree of dominance of T over C allele. Alike the CSN3, whose additive effect was significantly associated with SCS, a degree of dominance has been observed also for milk urea, where the heterozygous had significantly higher average values when compared to the opposite homozygous (Table 5).

Milk somatic cells consist of milk-secreting cells and immune cells. Regarding the CSN3, it is known that it derives from the fibrinogen by a duplication gene event [46] and that fibrinogen is one of the main mediators of inflammation acute phase [47]. Therefore, it is possible that the κ-casein kept part of the ancestor gene’s functions and plays an active role as indicator of SCS and mastitis. A further support to this statement derives from the function that the κ-casein glycomacropeptide (GMP) carries out in the modulation of immune response and as antibacterial and anti-inflammatory peptide [48,49,50]. In addition, recently in domestic cattle, the SNP rs43703017, located in the CSN3, was associated with an increase of SCS [51]. Regarding the CSN1S1, the association with SCS confirmed in buffalo the significant effect of this gene as interesting candidate for selection to improve resistance against mastitis as already indicated in dairy cows [52,53].

Looking at milk urea, no genes had significant substitution effect on this trait. The polymorphism on SCD seems to not affect any of the investigated milk phenotype for AM. Dominance positive effects are suggested (p<0.05) for SCS (CSN1S1) and milk urea (CSN3).

The use of the genotypic model substantially confirmed the results of allelic model with few differences in the significance level for LPL (dMY, dPP), αs1-CN (dPY) κ-CN (dPP) that only approached the significant threshold (p<0.10) but with a good approximation can be considered suggestive of a SNP-phenotype association as also confirmed by the proportion of variance explained by SNP effects for those trait-gene association (from 0.2% to 0.4%) (Table 5). Indeed, LPL polymorphism accounted for 0.3% and 0.2% of total variability dMY and dPP respectively. The polymorphism at CSN1S1 explained the 0.4% of total variance for dMY and dPY and SCS. Despite the reduced percentage of variance in absolute values (0.1% to 0.7% cumulatively across traits), this is not unusual when genetic association of single genes are analyzed.

It is worth noting that the random effect of buffalo cows and htd explained a large part of variance. In general, it appears that variance accounted for buffaloes are larger for SCS and urea (25%-57%) and smaller for milk yield and composition (8%-14%). With an opposite trend, htd largely explain intra-herd-test-day variability (26-37%) for dMY, dPY and dFY and less of milk contents, SCS and urea (7.5%-14.5%). In this context, the different environmental and management conditions among the eight farms might have not allowed for a better control of some sources of non-genetic variation. Therefore, such high level of variability observed in the present study may be ascribed to the relevant effect of environmental noise.

Representative examples of DIM classes least square means for dPP and dFP for the three LPL genotypes and dPP and SCS for CSN3 (Figure 1) are reported within lactation pattern of different genotypes.

5. Conclusions

The genetic improvement of dairy traits is among the main goals of the Italian Mediterranean river buffalo association of breeders. In the present study, we have extended to a larger population the investigation on four polymorphisms that, previously in a limited number of samples and often in a single farm, have been associated to dairy traits. Three out four SNP in CSN1S1, CSN3 and LPL were associated with at least one of analyzed traits (dMY, dPY, dFY, dPP, dFP, SCS and Urea) using both an allelic and a genotypic model. In particular, the CSN1S1 (AJ005430:c.578C>T) gave favorable associations with all yield traits (dMY, dPY, dFY) and SCS, whereas CSN3 positively associated with SCS and Urea. Favorable effect on dMY, dFP and dPP were observed for the LPL. Conversely, SCD did not show any association with milk traits. Overall, our results are important indications for the preselection of young buffalo bulls for the dairy traits, but they also highlight the importance of confirmation studies in larger populations to validate previous association studies for a more efficient setup of gene-assisted breeding programs.