The Italian buffalo dairy industry has increased in the last fifty years mainly owing to the growing economic importance of its main product, the mozzarella PDO (Protected Denomination of Origin – Reg. EC 510/2006) cheese. In the last decade the Italian buffalo stock has almost doubled (http://www.aia.it), with approximately 344 000 Mediterranean river buffaloes (Bubalus bubalis) farmed in Italy at present. Most of the herds, about 2500, are located in the south of the Country (http://agri.istat.it).
Such an increase of economic relevance has not being accompanied by an improvement of the production level of the animals. Average milk yield per buffalo cow is 2180±597 kg in 270 d of lactation (it was 2140 kg in the year 2000) with 8·45% fat and 4·59% protein (AIA, 2010). Although milk composition fits cheese processing requirements, the productive level is insufficient to satisfy the market demand and for meeting economic goals of farmers. Thus both management and breeding improvements are needed.
The Italian buffalo Herd book has 35 114 registered buffaloes. Most of them are involved in a dairy recording programme. However, the implementation of a conventional breeding programme for dairy traits based on progeny testing and EBV calculation is hampered by the poor efficiency of AI in this species. Difficult oestrus detection and variability of its length are among the main causes (Barile, Reference Barile2005; Drost, Reference Drost2007). In fact, natural mating is the most widely used technique of reproduction. Moreover, small values of estimated genetic parameters have been ascribed to inaccurate identification of true paternity (Rosati & Van Vleck, Reference Rosati and Van Vleck2002).
Marker-assisted selection may represent a possible option for designing a suitable breeding scheme for Italian river buffaloes. Gene polymorphisms significantly associated with milk production traits may provide useful indications for identifying selection candidates with high genetic merit. Actually, few studies have been carried out on candidate genes in buffaloes. One example is represented by the association between α-lactalbumin gene polymorphism and milk yield and composition (Dayal et al. Reference Dayal, Bhattacharya, Kaushik and Singh2009). More studies can be found on the molecular characterization of loci that may potentially affect dairy traits. Of particular interest is a polymorphism recently reported for the oxytocin (OXT) gene (Cosenza et al. Reference Cosenza, Pauciullo, Mancusi, Nicodemo, Di Palo, Zicarelli, Di Berardino and Ramunno2007). This locus, mapped on river buffalo chromosome14q22 (Iannuzzi et al. Reference Iannuzzi, Gallagher, Di Meo, Schlapfer, Perucatti, Amarante, Incarnato, Davis, Taylor and Womack2001), may be a possible candidate gene for improving milk yield and milkability, owing to the role of the oxytocin hormone in alveolar milk ejection and in milk flow rate. A successful milking requires a complete milk removal from both cisternal and alveolar compartments of the mammary gland. For a complete milk ejection, oxytocin must be released from the pituitary gland and transported to the udder where it acts on myoepithelial cells to promote contraction (Bruckmaier & Blum, Reference Bruckmaier and Blum1998). Milk ejection in response to suckling or milking is achieved via a classical neuro-endocrine process, known as the milk ejection reflex (Wakerley et al. Reference Wakerley, Clarke, Summerlee, Knobil and Neill1994). This process is of particular importance for buffaloes, in which alveolar milk represents about 95% of total milk owing to the absence or small size of the udder cistern (Thomas et al. Reference Thomas, Svennersten-Sjaunja, Bhosrekar and Bruckmaier2004; Caria et al. Reference Caria, Murgia and Pazzona2011).
The OXT gene is 912 bp long and it codes for 106 amino acids of the oxytocin complex. Three SNPs have been discovered in Italian river buffalo: two transitions in the promoter 5′flanking region (AM234538: g.28C>T and g.204A>G) and a non-synonymous transversion (g.1627G>T) at the 170th nucleotide of the second exon. The latter SNP is responsible for the Arg97→Leu amino acid substitution in the mature protein, which yields two alleles respectively named A (EMBL Acc. No AM234538) and B (EMBL Acc. No AM234539) (Cosenza et al. Reference Cosenza, Pauciullo, Mancusi, Nicodemo, Di Palo, Zicarelli, Di Berardino and Ramunno2007).
The aim of the present work was to test possible associations between the above mentioned SNP polymorphisms at the OXT locus and milk yield in the Italian Mediterranean river buffalo.
Materials and Methods
DNA sampling and genotyping
Biological samples were collected from 170 Italian river buffaloes belonging to an experimental herd, located in Salerno province (Southern Italy). Sampling was carried out in collaboration with the Italian National Association of Buffalo Breeders (ANASB). DNA was extracted from blood leucocytes using the procedure described by Gossens & Kan (Reference Gossens and Kan1981).
Briefly, 5 ml of frozen blood samples were washed twice with distilled water and 1.8% NaCl, protein digestion was carried out with 500 μl of proteinase K solution (2 mg/ml of proteinase K, 1% w/v SDS and 0·02 m-EDTA). Proteins were extracted using a phenol-chloroform method followed by DNA precipitation with cold isopropanol. The isolated DNA was then resuspended in 100 μl TE buffer pH 7·6 (10 mm-Tris–1 mm-EDTA).
DNA concentration and OD260/280 ratio of the samples were then measured with the Nanodrop ND-2000C Spectrophotometer (Thermo Scientific). The entire panel of 170 animals was genotyped for three SNPs in out-sourcing (http://kbioscience.co.uk) using a competitive allele specific PCR system (KASPar), a novel form of homogeneous allele-specific fluorescent genotyping system based on Fluorescent Resonance Energy Transfer technology (FRET) and the use of VIC/FAM fluorochromes.
Phenotypic data collection
A total of 41 980 test-day records belonging to 219 lactations of 163 buffalo cows of different parities (1–11) were used. Data were collected in the period January 2007–September 2009. Animals were milked twice a day with the Afifarm system (S.A.E. Afikim, Kibbutz Afikim, Israel). Milk yield and milk flow were recorded individually at each milking for the whole lactation through an Afi-Flo 2000 flow-meter. Only animals with complete genotyping for the three SNPs were retained for the statistical analysis. Average number of records per lactation was 192±97.
Statistical analysis
Allelic frequencies were calculated for each SNP locus. Measures of linkage disequilibrium (D′ and r 2) were estimated using Haploview software version 4·2 (http://www.broadinstitute.org/haploview/haploview). Haplotype structure was defined according to Gabriel et al. (Reference Gabriel, Schaffner, Nguyen, Moore, Roy, Blumenstiel, Higgins and Altshuler2002).
Associations between OXT polymorphisms for 163 buffaloes and milk yield and flow were investigated by fitting the following mixed linear model with the MIXED procedure of SAS 9.2 (2008):
where y ijklmno=milk yield or milk flow in test o; μ=overall mean; Yeari=fixed effect of the ith year of production (3); Parj=fixed effect of the jth parity (6 classes: 1 to 5, >5); Seak=fixed effect of the kth season (4); OXTl=fixed effect of the lth OXT genotype; DIMm=fixed effect of the mth stage of lactation (30 intervals of 10 d each); c n=random effects of individual buffalo cow nested within oxytocin genotype e ijklmno=random residual.
The fixed effect of SNP genotype fits mean gene effect across the whole lactation. The days in milk (DIM) factor nested within OXT genotype was included in the model to estimate lactation curves of the different genotypes (Stanton et al. Reference Stanton, Jones, Everett and Kachman1992). (Co)variance matrices of random effects of buffalo cow and residual were assumed to be diagonal, Iσ2c and Iσ2e, respectively. They allow for the REML estimation of variance components associated to individual buffalo cow (σ2c) and residual (σ2e). Statistical significance of the SNP effect was tested against variance of buffalo cow nested within SNP genotype according to repeated measures design theory (Littell et al. Reference Littell, Henry and Ammerman1998). Pairwise comparisons among different levels of fixed effects included in the model were performed using a Bonferroni adjusted test.
The average gene substitution effect (α) was calculated using a mixed linear model with the same structure as (Eq. 1) but with the gene effect treated as a covariable, represented by the number of G alleles at the OXT locus (0, 1, 2), and an interaction between alleles at the SNP locus to account for possible dominance effects (Banos et al. Reference Banos, Woolliams, Woodward, Forbes and Coffey2008).
Finally, in order to estimate the contribution of the oxytocin locus to the variance of the trait, a mixed model having the same structure of (Eq. 1) but with the OXT genotype treated as random was run. Thus a variance component associated with the OXT locus (σ2OXT) was estimated. Contributions of OXT locus (r 2OXT) and buffalo cow (r 2c) to the total phenotypic variance of the trait considered were calculated as:
Results
Genotype and allele frequencies together with distribution of buffalo cows, lactations and tests across different SNP genotypes are reported in Table 1. On average, the allele call rate was 98·8%. In all the three loci the major allele had a relative frequency of about 0·76 and χ 2 values showed that there was no evidence for departure from the Hardy-Weinberg equilibrium (P⩽0·05).
Genotyping data, allele frequencies, Hardy–Weinberg equilibrium (P⩽0·05), and absolute frequencies of buffalo cows, lactations and tests across genotypes of the three SNPs at the OXT locus
A complete linkage disequilibrium was observed. Average D′ and r 2 values were 1 and 0·97 respectively. Only two haplotypes were found: CAG and TGT with frequencies of 0·76 and 0·24 respectively, i.e. equal to allelic ones. According to Gabriel et al. (Reference Gabriel, Schaffner, Nguyen, Moore, Roy, Blumenstiel, Higgins and Altshuler2002), loci can be considered as being in the same haplotype block. Thus only one SNP (g.1627G>T) was considered in running model (Eq. 1).
Both milk yield and flow were significantly affected by all fixed factors included in model (Eq. 1) (Table 2), except for OXT genotype for milk flow (Table 3). Milk yield tended to increase from first to latter parities, reaching the maximum at the fourth calving. This parity order showed also the highest value for milk flow. Actually the two variables are highly correlated (r=0·97). Largest daily milk yields were observed for buffaloes calving in autumn whereas lowest values were for summer calvings (Table 2). Moreover milk flow showed the highest values in the autumn calvings.
Least squares means of milk yield (kg/d) and milking speed (ml/s) for the different levels of parity and calving season estimated with model (Eq. 1)
A,B Means within columns without a common superscript differ (P<0·01)
a,b Means within columns without a common superscript differ (P<0·05)
Least squares means of milk yield (kg/d) and milking speed (ml/s) for the genotypes at the locus g.1627G>T of river buffalo oxytocin gene estimated with model (Eq. 1)
a,b Means within columns without a common superscript differ (P=0·029)
Contribution of the individual buffalo cow effect to the total phenotypic variance, also named repeatability (Stanton et al. Reference Stanton, Jones, Everett and Kachman1992), was about 0·57 for milk yield and 0·10 for milk flow, respectively (Table 4).
Substitution effect of a guanine for a thymine at the locus g.1627G>T of the oxytocin gene (mean±se) and contribution of the OXT locus at the total phenotypic variance
† α: Substitution effect; d: dominance effect; σ2: variance components associated to the genotype (OXT); to the individual buffalo cow (c), to residuals (e); r 2: contributions of genotype (OXT) and of individual buffalo cow (c) to the total phenotypic variance
The OXT genotype was significantly associated with milk yield (P=0·029) (Table 3). In particular, TT buffaloes showed an average daily milk yield approximately 1·7 kg higher than GT buffaloes (Table 3). Such a relevant difference, representing about 23% more milk/d, started at approximately 60 d of lactation and tended to remain constant throughout the whole lactation as it can be seen from the estimated lactation curves of the three genotypes reported in Figure 1. Compared with homozygous genotypes, GT buffaloes had a slightly lower average daily yield. A large dominance effect, statistically different from zero, was estimated. It is double the average substitution effect (Table 4). The contribution of OXT genotype (r 2OXT) to the total phenotypic variance in milk yield was equal to 0·01 (Table 4).
Lactation curves of different OXT genotypes for milk yield (kg/d). The GG genotype comprised 94 buffalo cows and 123 lactations, the GT 60 buffalo cows and 82 lactations, the TT 9 buffalo cows and 14 lactations, respectively.
The TT genotype showed highest values also for milk flow (20·8 and 14·7% more milk flow compared with GT and GG genotype respectively), even though the estimated difference did not reach a level of significance (P=0·07) (Table 3). However, as for milk yield, also the milk flow curve was higher across the whole lactation for TT buffaloes (Fig. 2) compared with other genotypes. Also for this trait, an over-dominance was evident (Table 3). The contribution of OXT genotype (r 2OXT) to the total phenotypic variance in milk flow is lower than in the case of milk. This trend can be observed also for the individual contribution of the buffalo cow effect (Table 4).
Lactation curves of different OXT genotypes for milking flow (ml/s). The GG genotype comprised 94 buffalo cows and 123 lactations, the GT 60 buffalo cows and 82 lactations, the TT 9 buffalo cows and 14 lactations, respectively.
Discussion
In this work, an association between three SNPs at the OXT locus and daily milk yield of Italian river buffalo is reported. Although the results refer to a single herd and should be tested population-wide, it is of great importance as one of the first indications of an association between a trait of economic importance and a candidate locus in this species. Moreover, it is among the first reports of an association between the oxytocin locus and production traits in ruminants. Therefore, it may offer useful indications for marker-assisted selection not only in buffaloes, but also for starting investigations in other dairy species.
Environmental factors affecting milk traits
Buffalo daily milk yield and flow were found to be significantly affected by parity, calving season and month of milking (Table 2). As expected, primiparous buffaloes showed a lower milk yield. The observed higher milk yield for third and fourth parities is in agreement with previous reports on Italian (Catillo et al. Reference Catillo, Macciotta, Carretta and Cappio-Borlino2002) and Murrah (Dang et al. Reference Dang, Mukherjee, Kapila, Mohanty, Kapila and Prasad2010) buffaloes. Moreover the milk flow tended to increase from first to later parities, even though with a less definite pattern, in agreement with results reported for dairy cattle (Pfeilsticker et al. Reference Pfeilsticker, Bruckmaier and Blum1996; Sandrucci et al. Reference Sandrucci, Tamburini, Bava and Zucali2007; Guler et al. Reference Guler, Yanar, Aydin, Bayram, Dogru and Kopuzlu2009). An effect of parity on milk flow and total milk volume has been reported also for the goat, where second and third lactation animals show higher values (Peris et al. Reference Peris, Such and Caja1996), whereas conflicting results have been observed for dairy sheep (Dell'Aquila et al. Reference Dell'Aquila, Pilla, Catillo, Scardella and Taibi1993; Sevi et al. Reference Sevi, Taibi, Albenzio, Muscio and Dell'Aquila2000).
The relevant influence of seasonal effects on productive performances of buffaloes has been confirmed in the present work. The lowest average yield for calving in summer can be explained mainly by the depressive effect of high temperatures at the beginning of the lactation. Similar figures were obtained in Mediterranean and Pakistan Nili-Ravi breeds (Catillo et al. Reference Catillo, Macciotta, Carretta and Cappio-Borlino2002; Afzal et al. Reference Afzal, Anwar and Mirza2007). The effect of calving season observed for milk flow, with highest flows for buffaloes that calve in autumn, is different from that reported for dairy cattle (Firk et al. Reference Firk, Stamer, Junge and Krieter2002), and is probably due to the physiological and management differences between both species.
Association between OXT genotype and milk traits
The three SNPs at the OXT gene investigated in this study were found to be in the same haplotypic block. This was an expected result, given their close distance (only 1136 bp separate the first and the last polymorphic site). The transversion g.1627G>T at the 170th nucleotide of the second exon responsible for the Arg97→Leu amino acid change was chosen for the association study. This SNP is in the same linkage phase with the g.204A>G polymorphism in the promoter region of OXT gene, which is responsible for the presence/absence of a potential binding site for CCAAT/enhancer binding protein-alpha (C/EBP alpha) (by TRANSFACT 7.0 database) (Cosenza et al. Reference Cosenza, Pauciullo, Mancusi, Nicodemo, Di Palo, Zicarelli, Di Berardino and Ramunno2007). The C/EBP alpha cooperates with other transcription factors, including pituitary-specific homeodomain protein (Pit-1), in the regulation of lactotropes gene transcription in the anterior pituitary. An effect of promoter polymorphism on milk production traits has been reported in cattle for the DGAT1 gene by Fürbass et al. (Reference Fürbass, Winter, Fries and Kühn2006). Although a general action of enhancement of transcriptional activity might offer a possible explanation of the action of such OXT polymorphism on milk yield, further investigations are necessary.
The superiority in milk yield showed by TT buffaloes over the heterozygous genotype (more than 1·7 kg/d) is relevant compared with the relatively low level of production of this species. An interesting feature of such a difference is that it remains quite constant throughout the whole lactation. This feature has been observed also for other candidate genes found to affect dairy traits in cattle such as DGAT1 (Strucken et al. Reference Strucken, de Koning, Rahmatalla and Brockmann2011), where the effect can be observed after 40 DIM, and SCD (Macciotta et al. Reference Macciotta, Mele, Conte, Serra, Cassandro, Dal Zotto, Cappio Borlino, Pagnacco and Secchiari2008). The remarkable magnitude of the OXT effect is confirmed by its contribution to the total phenotypic variance of milk yield. The observed value is larger than the one reported for the DGAT1 (Grisart et al. Reference Grisart, Coppieters, Farnir, Karim, Ford, Berzi, Cambisano, Mni, Reid, Simon, Spelman, Georges and Snell2002) but less than half of the contribution reported for the ABCG2 (Cohen-Zinder et al. Reference Cohen-Zinder, Seroussi, Larkin, Loor, Everts-van der Wind, Lee, Drackley, Band, Hernandez, Shani, Lewin, Weller and Ron2005) genes on milk yield in dairy cattle.
The marked unbalance of genotype class data, with the best genotype represented by less than 10 buffalo cows demands great care in making inferences from the observed results. Actually, outliers may have a relevant impact on the estimated effect of the gene. Also the location of the investigated animals in a single herd should be mentioned. However, some elements have to be considered as possible counterbalances to the above mentioned drawbacks. Firstly, the large number of observations available per buffalo cow (about 200) that might have positively affected the accuracy of estimated differences. Secondly, the strictness of the statistical tests used. The approach used for testing the effect of the OXT prevents spurious statistical significances that may arise when the residual variance is used in denominator of the F test. Moreover, the Bonferroni correction for the comparisons between the different genotype classes is very conservative. Finally, the common environmental and management conditions shared by the buffaloes might have allowed for a better control of some sources of non-genetic variation.
A large effect of dominance on milk yield was observed in the present work. In dairy cattle, dominance effects for milk production traits have been reported (Weller et al. Reference Weller, Golik, Seroussi, Ezra and Ron2003; Kuehn et al. Reference Kuehn, Edel, Weikard and Thaller2007). Often dominant allele effects are not detected or considered not relevant, being numerically much lower than the additive effects, even though they might have an impact on allele substitution effects in the population.
Results for milk flow are of a lower magnitude as far the effect of the gene is concerned. Its contribution to the total phenotypic variance in milk flow is markedly lower than in the case of milk. This trend can be observed also for the individual contribution of the buffalo cow effect. Such a result is quite unexpected, considering the physiological role of oxytocin and the correlation between milk yield and ejection. However, it should be remembered that previous studies carried out in buffalo and cow highlight that it is the increase of the oxytocin level over a threshold that affects milk ejection and not its absolute concentration (Schams et al. Reference Schams, Mayer, Prokopp and Worstorff1984; Thomas et al. Reference Thomas, Bruckmaier, Östensson and Svennersten-Sjaunja2005a). Thus a possible explanation may be that OXT genotype influences the total content of the hormone, but not its increase above the basal level. In any case, also the greater variability of milk ejection compared with milk yield (coefficient of variability 0·72 and 0·43, respectively) may be responsible for the lack of the statistical significance of the OXT genotype for this trait.
The estimated correlation among test-day milk yield within buffalo cow agrees with previous reports obtained on Italian river buffaloes (Catillo et al. Reference Catillo, Macciotta, Carretta and Cappio-Borlino2002), whereas no information is available yet far for milk flow in such species. Therefore, the value for milk flow can be compared with results obtained in other species. In cattle, the repeatability of milk flow measurements has been estimated as 0·54 and 0·59 for Brown Swiss and Simmental breed, respectively (Ilahi & Kadarmideen Reference Ilahi and Kadarmideen2004), whereas in dairy sheep it was 0·50 (Marie-Etancelin et al. Reference Marie-Etancelin, Casu, Aurel, Barillet, Carta, Deiana, Jacquin, Pailler, Porte, Tolu, Gabiña and Sanna2002). The lower value observed in the present work may be ascribed to the relevant effect of environmental noise. Buffaloes, in fact, are very sensitive even to small changes in the milking routine, and disturbed milk ejection occurs frequently (Thomas et al. Reference Thomas, Bruckmaier, Östensson and Svennersten-Sjaunja2005a, Reference Thomas, Nordstrom, Svennersten-Sjaunja and Wiktorssonb). Accordingly, measurements within buffalo cow were not strongly related.
Conclusions
In recent years, a great deal of work has been done on Italian buffalo farms to improve recording, health, feeding and livestock systems. On the other hand, little has been done in terms of genetic improvement in this species.
In the present study, novel genotyping information and haplotype structure of a key gene involved in milk production were reported. A significant genetic association between three SNPs at the oxytocin locus and daily milk yield of Italian river buffalo was found. The genotype with selective advantage showed an average daily milk yield approximately 1·7 kg higher. Such a result, if confirmed population-wide, may be of great economic interest for the buffalo dairy industry. Increases in average milk yield and, consequently, in mozzarella PDO production may be expected.
The same polymorphisms showed a less defined effect on milk flow. Results suggest a need for further investigations on a wide number of herds. Milk production and milk flow were significantly affected also by parity, and calving season, confirming the seasonality of the productive performances of such species.
To our knowledge, this is one of the first indications of genetic association between a trait of economic importance in buffalo. Furthermore, it is among the first reports of an association between the oxytocin locus and production traits in ruminants. Therefore, it may offer useful indications for further analyses and application of marker-assisted selection programmes in other dairy ruminant species.
This work was financially supported by the Italian Ministry for Agriculture and Forestry Policy – MiPAAF (SELMOL project).
