Animals and phenotypic records

Blood samples were collected from 465 Italian Mediterranean river buffalo belonging to four experimental farms in Southern Italy. The total milk production records used in the present study were obtained from Italian Buffalo Breeders Association (ANASB), which is responsible for recording the production of buffalo milk. A total of 1214 lactation records including milk yield (total lactation yield and also peak yield), protein content and fat content were collected. All the milk production records were adjusted to 270 d (Baldi et al. Reference Baldi, Laureano, Gordo, Bignardi, Borquis, de Albuquerque and Tonhati2011). The Ethical Animal Care and Use Committee of Federico II University of Naples (Italy) approved the experimental design and animal treatment protocols.

Buffalo genomic DNA was extracted from whole-blood samples using the QIAamp DNA Blood kit (Qiagen, Milano, Italy). The purity and concentration of DNA was measured by spectrophotometric, NanoDrop 2000 analyser (Thermo Fisher Scientific, Wilmington).

Single nucleotide polymorphism (SNP) identification and genotyping

The SNP primers were designed based on GenBank sequence (GenBank accession number: 102412882) for the PRL gene by primer premier5·0 software (Premier Biosoft, Palo Alto, CA) and synthesised by Quintara Biosciences, Wuhan. The identification of the SNP was performed by sequencing.

A high resolution melting (HRM) method was used to genotype all samples with the Light Cycler 480 System (Roche, no. 04909631001). In order to obtain great accuracy of genotyping in HRM analysis, new primers (Supplementary Table S1) were designed to amplify the small fragments (100–400 bp) including the SNPs. The reaction mixture included genomic DNA (20 ng), HRM Master Mix (2 × conc.) 10 µl, Primer mix (20 × conc.) 1 µl, MgCl₂ (25 mm) 2 µl and water (PCR-grade) to a total volume of 20 µl. Cycling conditions included an initial 10 min at 95 °C, followed by 50 cycles at 95 °C for 10 s, annealing temperatures (Supplementary Table S1) for 15 s, and 72 °C for 20 s, a single cycle of 95 °C for 1 s, 72 °C for 90 s, and melting tempering increasing from 70 to 90 °C at 0·1 °C/s. Gene Scanning Software was applied for analysis the HRM data. The different genotypes were determined based on melting curves of control samples.

Statistical analysis

Allelic frequencies, genotypic frequencies, Linkage disequilibrium, haplotypes and Hardy-Weinberg equilibrium testing were calculated for each SNP locus by the SHEsis online platform. Association between PRL polymorphisms for 465 buffaloes and milk production traits were investigated by fitting the following mixed linear model with the PROC MIXED procedure of SAS9·2:

(1) $$y_{ij} = \mu + Par + Sea + PRL + C_i (PRL) + \varepsilon _{ij} $$

where y _ij is milk production traits; μ is the overall mean; Par is the effect of the parity (five classes: 1–4, >4); Sea is the fixed effect of the season (four classes); C _i is the random effect of individual buffalo cow nested within the PRL genotype; and ε _ij is the random residual.

The fixed effect of the SNP genotype fits the mean gene effect across the whole lactation. In order to estimate the contribution of the SNP locus to the variance of the traits, a mixed model having the same structure of Eq. (1) but with the PRL genotype treated as random was run. This model, including a variance component associated with the PRL locus (δ² _PRL ) was estimated. Contributions of the PRL locus (r ² _PRL ) and the individual buffalo cow (r ² _c) to the total phenotypic variance were calculated as:

$$r_{PRL}^{\rm 2} = \delta _{PRL}^2 /(\delta _{PRL}^2 + \delta _{\rm c}^2 + \delta _{\rm e}^2 )$$

and

$$r_{\rm c}^{\rm 2} = \delta _{\rm c}^2 /(\delta _{\rm c}^2 + \delta _{\rm e}^2 + \delta _{PRL}^2 )$$

where δ² _PRL is a variance component associated with the PRL locus; δ² _c is the variance from individual buffalo cow; δ² _e is the error variance.

SNP identification and genotyping

The results of sequencing identified four variants in intron1 (2609C > T), exon2 (2836C > T), intron2 (4621A > G) and exon5 (8434C > T) (Supplementary Fig. S1). Here, only exon2 (2836C > T) was a nonsynonymous switch, resulting in an amino acid change from arginine (Arg) to cysteine (Cys) at position 12. The HRM analysis showed a difference between the temperatures of melting (T _m) for the heterozygote, homozygous mutation and wild type as different genotypes due to different bonds (Supplementary Fig. S2). Also, direct sequencing showed the genotypes are in agreement with the HRM analysis. The genotypic and allelic frequencies are given in Supplementary Table S2. The genotypic frequencies of the polymorphisms were observed to be in Hardy–Weinberg equilibrium (P > 0·05), indicating they were influenced less by the selection pressure. High linkage disequilibrium was observed. Average D′ and r ² values were 0·767 and 0·408 respectively. Three haplotypes were found: CACT, CATT and TGCC with frequencies of 0·04, 0·87 and 0·03 respectively.

The well-known A–G mutation in the codon for amino acid 103 in exon3 of PRL which produces a polymorphic RsaI site (Lewin et al. Reference Lewin, Schmitt, Hubert, van Eijk and Arnheim1992) in dairy cattle, was not found in the present study, Similar results have been reported in the Khuzestan-Iran buffalo (Tabar et al. Reference Tabar, Fayazi, Roshanfekr, Mirzadeh and Sadr2012) and the Murrah Buffalo (Biradar et al. Reference Biradar, Unaune, Dodamani, Mhatre, Londhe, Pawar, Sawane and Umrikar2014), indicating the monomorphic nature of the locus in buffalo species. However, the amount of milk and fat yield of the GG genotypes was higher compared to the AA genotype in PRL-exon3 and the AA genotype was negatively correlated to the milk fat yield in black-and-white cattle (Khatami et al. Reference Khatami, Lazebny, Maksimenko and Sulimova2005). This difference between the buffalo and dairy cattle may be one of the reasons for the low milk yield and high milk fat and milk protein content in buffalo.

Association analysis

This is the first report to correlate disruption in the function of the PRL gene with changes in overall performance in terms of milk yield, fat and protein in buffalo. A total of 15 variations of PRL gene have been identified in Pakistan buffalo (Nadeem & Maryam, Reference Nadeem and Maryam2016), but no correlation analysis with milk yield trait was reported. In the present study, SNPs of intron2 and exon5 were not associated with milk production traits in the population. However, the SNPs of intron1 and exon2 showed significant associations with milk production traits, and were considered in a running model (Eq. (1)) Moreover, milk traits were affected significantly by all fixed factors included in the model (Eq. 1; Table 1).

Table 1. Least-squares mean of 270 d milk yield, protein content, fat content and peak milk yield for different levels of parity, calving season and genotypes at the locus intron1 (left) and exon2 (right) of the river buffalo PRL gene estimated with model Eq. (1)

^a,bMeans within a row with different superscripts differ significantly (P < 0·05)

Our findings showed that intron1 polymorphism (2609C > T) was associated significantly (P < 0·05) with 270 d milk yield, total protein yield and milk peak (Table 1). The milk yield tended to increase from the first to latter parities, reaching a maximum at the fourth calving. The parity order showed the highest value for milk protein content and milk yield. Buffalo which calved in autumn had the highest milk yield. By contrast, the smallest milk yield was observed for buffalo calving in summer, which showed the lowest values of all other milk production traits. In particular, the milk production of buffalo with TT genotype was 49 kg less than other genotypes of buffalo in 270 d. In addition, buffalo with TT genotype showed the lowest value of total protein yield and peak milk yield (P < 0·05) when compared with buffalo with CT or CC genotype. The contribution of the individual buffalo effect on milk yield, protein content and peak milk yield was about 0·21, 0·02 and 0·06, respectively. The average effect of the intron1 genotype to the total phenotypic variance on milk production traits was 0·09 (Table 2).

Table 2. The polymorphism effect at the locus intron1 of the PRL gene (mean ± se) and contribution to the total phenotypic variance

The transition g.2836C > T at 12^th nucleotide of the second exon responsible for arginine to cysteine change showed an important effect on milk fat content. The milk fat content tended to increase from first to later parities, reaching the maximum at the fifth calving and the buffalo which calved in winter had the highest value of milk fat content. In addition, g.2836C > T polymorphism was associated (P < 0·05) with the largest total fat yield, the buffalo with TT genotype generally had better value of total fat yield than CT genotype groups. No CC group was observed in the studied population (Table 1).