The buffalo is a relatively young species, domesticated approximately 5000 years ago. The little exchange of buffaloes among countries allowed a good differentiation among the known breeds, thus nowadays each population has maintained its own phenotypic traits and performances (Borghese & Mazzi, Reference Borghese and Mazzi2005). A long period of isolation and the lack of crossbreeding allowed an evident morpho-functional differentiation of the Mediterranean type, which has been recognized as ‘Mediterranean Italian’ breed since 2000 (D.M. 201992 – 5/7/2000).
The Italian buffalo population has increased of 53·5% in the decade 1999–2009 (http://www.aia.it). The reason for this trend lies in the growing interest of farmers in buffalo milk production (not subject to EU milk quotas), and the increasing market demand for mozzarella PDO (Protected Denomination of Origin—Reg. EC 510/2006) cheese. Such an increase of economic importance and size of the dairy Buffalo industry was followed by a limited 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 and 4·59% of fat and protein contents respectively (AIA, 2010).
Gene polymorphisms significantly associated with milk production traits may provide useful indications for identifying selection candidates with high genetic merit. In this regard, a single nucleotide polymorphism (SNP) (EMBL acc. no. FM876222: g.133A > C) recently reported in the promoter region of the river buffalo stearoyl-CoA desaturase (SCD) gene (Pauciullo et al. Reference Pauciullo, Cosenza, D'Avino, Colimoro, Nicodemo, Coletta, Feligini, Marchitelli, Di Berardino and Ramunno2010) is particularly interesting.
Stearoyl-CoA desaturase is a microsomal enzyme which plays a key role in fatty acid metabolism. It is also known as delta-9-desaturase because it catalyses the introduction of the first cis-double bond in the Δ9 position in a large spectrum of fatty acyl-CoA substrates (Ntambi, Reference Ntambi1999). The stearoyl-CoA desaturase (SCD) locus has been suggested as a candidate gene affecting milk fatty acid (FA) profile (Gautier et al. Reference Gautier, Barcelona, Fritz, Grohs, Druet, Boichard, Eggen and Meuwissen2006). In river buffalo, the g.133A > C SNP creates a new consensus site for the SP1 binding site and a preliminary association study with the milk FA content showed a significant effect on the desaturation index (Pauciullo et al. Reference Pauciullo, Cosenza, D'Avino, Colimoro, Nicodemo, Coletta, Feligini, Marchitelli, Di Berardino and Ramunno2010). In sheep, no significant associations have been reported (García-Fernández et al. Reference García-Fernández, Gutiérrez-Gil, García-Gámez, Sánchez and Arranz2010), even though this gene is a strong positional and functional candidate for a quantitative trait locus (QTL) detected on chromosome 22 for the conjugated linoleic acid (CLA)/vaccenic ratio in milk (Carta et al. Reference Carta, Casu, Usai, Addis, Fiori, Fraghı, Miari, Mura, Piredda, Schibler, Sechi, Elsen and Barillet2008). Associations between polymorphisms in SCD gene and the FA composition of milk have been reported in cattle (Mele et al. Reference Mele, Conte, Castiglioni, Chessa, Macciotta, Serra, Buccioni, Pagnacco and Secchiari2007; Moioli et al. Reference Moioli, Contarini, Avalli, Catillo, Orru, De Matteis, Masoero and Napolitano2007); and, recently a polymorphism at bovine SCD locus has also been associated with daily milk yield with a constant effect of the genotype across lactation stages (Macciotta et al. Reference Macciotta, Mele, Conte, Serra, Cassandro, Dal Zotto, Cappio, Pagnacco and Secchiari2008). Therefore, SCD locus has been proposed in gene-assisted selection programmes for the improvement of milk production traits in dairy cattle (Macciotta et al. Reference Macciotta, Mele, Conte, Serra, Cassandro, Dal Zotto, Cappio, Pagnacco and Secchiari2008).
The aim of this study was to investigate the occurrence of the above mentioned polymorphism at the SCD locus in the Italian Mediterranean river buffalo and to test possible associations with milk yield.
Materials and Methods
Sample collection and nucleic acid isolation
A total of 322 buffaloes from different commercial herds located in two provinces (Salerno and Caserta) of the Campania region (Southern Italy) were considered in the study. Most of the farmed buffaloes of the country are in this region of Italy.
Biological samples of these animals were used to extract DNA from blood, using the procedure described by Gossens & Kan (Reference Gossens and Kan1981).
DNA concentration and OD260/280 ratio of the samples were measured with the Nanodrop ND-2000C Spectrophotometer (Thermo Scientific).
PCR conditions and genotyping by TaqI PCR-RFLP
PCR reaction mixture and thermal conditions for the amplification of the DNA fragment spanning from −593 bp of 5′ flanking region to +159 bp of exon 1 (752 bp) of the river buffalo SCD gene were accomplished by using the following primer: SCD AF (5′-GAAACTTCCCCAGTGCC-3′) and SCD AR (5′-CAAGTGGGCCGGCATC-3′), according to Pauciullo et al. (Reference Pauciullo, Cosenza, D'Avino, Colimoro, Nicodemo, Coletta, Feligini, Marchitelli, Di Berardino and Ramunno2010). Product specificity was confirmed by ethidium-bromide-stained 1·5% agarose gel electrophoresis.
The entire panel of 322 animals was genotyped for the g.133A > C SNP using a PCR-RFLP method. Digestion of 17 μl of each PCR amplification was accomplished with 10 U of TaqI endonuclease (T↓CGA) (Promega) for 5 h at 65 °C. The digestion products were analysed directly by electrophoresis in 1·5% agarose gel in 1X TBE buffer and stained with ethidium bromide.
Phenotypic data collection
For a subsample of 169 buffalo cows an association study between the SNP polymorphism and milk yield was carried out. A total of 43 510 test-day (TD) records for milk yield measured daily with an automatized milk recording system on 226 lactations of 169 Italian water buffalo cows of different parities (1–7) were used. Data were collected in the period January 2007–September 2009 in a herd located in the province of Salerno. Animals were milked twice a day with the Afifarm system (S.A.E. Afikim, Kibbutz Afikim, Israel).
Bioinformatic and statistical analysis
The allele frequency and Hardy-Weinberg equilibrium (χ2 test) were calculated for the larger sample of 322 buffaloes. Homology searches, comparison among sequences, and multiple alignments were accomplished using DNAsis-Pro (Hitachi), whereas the putative transcription factor binding sites were searched by Transfact® 7.0 software.
Associations between SCD polymorphism for 169 buffaloes and milk yield was investigated with the following mixed linear model (SAS Institute, Cary NC, USA):
$$y_{ijlmn} = {\rm YS}_i + {\rm Par}_j + {\rm SCD}_l + {\rm DIM}_m + {\rm DIM}_m ({\rm SCD}_l ) + c_n ({\rm SCD}_l ) + {\rm} e_{ijlmn} $$where y ijlmn = test-day record of milk yield; YS = fixed effect of the ith year per season of production (8 levels: autumn 2007 to summer 2009); Par = fixed effect of the jth parity (5 levels: 1–4, >5); SCD = fixed effect of the lth SCD genotype (3 levels: AA, AC, CC); DIM = fixed effect of the mth stage of lactation (30 levels of 10 d each); c = random effects of individual cow (169 levels), nested within SCD genotype; e ijlmn = random residual.
The fixed effect of SNP genotype fits the average effect across the whole lactation. The DIM (days in milk) factor nested within SCD 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 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 cow (σ2c) and residual (σ2e). Statistical significance of the SNP effect was tested against variance of cow nested within SNP genotype (Littell et al. Reference Littell, Henry and Ammerman1998). Pairwise comparisons among different levels of fixed effects included in model were performed using a Bonferroni adjusted test.
The average gene substitution effect (α) was calculated using a mixed linear model with the same structure of (1) but with the gene effect treated as a covariable, represented by the number of C alleles at the SCD 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 SCD locus to the variance of the trait, a mixed model having the same structure of (1) but with the SCD genotype treated as random was run. Thus a variance component associated to the SCD locus (σ2SCD) was estimated. Contributions of SCD locus (r 2SCD) and cow (r 2c) to the total phenotypic variance of the trait considered were calculated as:
$$r^2 _{{\rm SCD}} = \displaystyle{{\sigma ^2 _{{\rm SCD}}} \over {\sigma ^2 _{{\rm SCD}} + \sigma ^2 _c + \sigma ^2 _e}} \quad{\rm and}\quad r^2 _C = \displaystyle{{\sigma ^2 _C} \over {\sigma ^2 _c + \sigma ^2 _e + \sigma ^2 _{{\rm SCD}}}} $$Results and Discussion
Characterization of river buffalo SCD promoter
A total of 593 bp of the river buffalo SCD promoter (EMBL acc. no. FM876222) upstream the first nucleotide of the first exon were analysed by using Transfact® 7.0 software and characterized in order to investigate the putative transcription factor binding sites that could regulate the gene expression (Fig. 1).
Key transcription factors binding sites found in the 5′ promoter region of the river buffalo SCD gene. The conserved polyunsaturated fatty acid (PUFA) response region including the sterol response element (SREBP), CCAAT-box (C/EBP), nuclear factor (NF)-1 and stimulator protein 1 (SP1) binding site are shown. TATA motifs and peroxisome proliferator activated receptor-γ (PPAR-γ) are also shown proximal to the transcription start site. SNP g.133A > C is indicated with M nucleotide according to international nomenclature
The buffalo SCD promoter contains two TATA box located, with reference to the first nucleotide of the first exon, at nucleotides −51/ − 48 and −102/ − 99 and at least four SP1 consensus sequence (−473/ − 468, −459/448, −58/64, −18/ − 6). The high GC content (64·9%) is a typical feature of housekeeping genes promoters (Zhu et al. Reference Zhu, He, Hu and Yu2008) and the presence of these last sites could contribute to the constitutive expression of the SCD in various tissues.
Several CCAAT/enhancer binding protein (C/EBP) were identified. This family of transcription factors are key regulators of adipogenesis and lipid metabolism, playing a fundamental role in expression of adipocyte genes as ADIPOQ, DGAT1, LPL, CD36 (Olofsson et al. Reference Olofsson, Orho-Melander, William-Olsson, Sjöholm, Sjöström, Groop, Carlsson, Carlsson and Olsson2008). In river buffalo SCD promoter, most of these C/EBP consensus sequences were found in a very closed DNA fragment of about 130 bp (nucleotide −382/ − 250), suggesting that this region could have an essential function in the gene expression. To support these findings, a multiple sequence alignment of the river buffalo, pig, human and mouse SCD − 1/SCD − 2 promoters was achieved (data not shown). Such DNA region was found to be highly conserved and among the found binding sites for transcription factors, those regulating lipid metabolism were: one sterol regulatory element binding protein (SREBP) at position −382/ − 372, which is considered a key components of trans-10, cis-12 C18:2 conjugated linoleic acid (CLA) regulation of bovine milk fat synthesis (Harvatine & Bauman, Reference Harvatine and Bauman2006); and two nuclear factor (NF)-1 binding sites (−362/ − 358 and −343/ − 339), which are functional enhancers with specificity for adipose cells and involved in adipocyte-specific gene expression (Graves et al. Reference Graves, Tontonoz, Ross and Spiegelman1991).
The analysis of the remaining portion of 5′ flanking region showed other regulatory elements: one sterol regulatory element (SREBP) at the nucleotides −71/ − 68, one NF-1 at the position −112/ − 108 and two peroxisome proliferator activated receptor-γ (PPAR-γ) at the positions −34/ − 29 and 10/15. PPAR-γ is often referred as the ‘master regulator’ of adipogenesis because it participates in the transcriptional activation of numerous adipogenic and lipogenic genes (Paton & Ntambi, Reference Paton and Ntambi2009). Recently, a gene network analysis in bovine mammary tissue showed that expression of PPAR-γ and its putative target genes was up-regulated during lactation, suggesting a role for this nuclear receptor in the regulation of milk fat synthesis (Bionaz & Loor, Reference Bionaz and Loor2008).
Genotyping
The transvertion g.133A > C creates a restriction site for the endonuclease TaqI, thus a PCR-RFLP protocol was set up for the quick genotyping of the samples. Digestion of the PCR product (752 bp) allows the identification of both alleles (Fig. 2).
Genotyping of SCD g.133A > C SNP by Taq I (T↓CGA) PCR-RFLP. Line 1: AA homozygous samples; line 3: CC homozygous samples; line 2 heterozygous samples. In the lines 2 and 3 is not visible the band 132 bp long. Line M is 2-log DNA ladder (0·1–10 kb; New England Biolabs)
The restriction pattern is characterized by one undigested fragment of 752 bp for the AA homozygous samples, whereas the same amplicon is restricted into two fragments of 132 and 620 bp in the presence of cytosine at the homozygous status. The restriction pattern of the heterozygous samples shows 3 fragments. The frequency of cytosine in the sample of 322 buffaloes was 0·16 (Table 1). This value might open the possibility for a rapid directional selection in favour of the C allele, which was found to be associated with a higher content of unsaturated FAs in milk (Pauciullo et al. Reference Pauciullo, Cosenza, D'Avino, Colimoro, Nicodemo, Coletta, Feligini, Marchitelli, Di Berardino and Ramunno2010). Hardy-Weinberg disequilibrium was detected for the genotype distribution, but this event could be the result of several causes. In particular, the restricted number of genotyped individuals and the mating system. In the last case, even though the mating technique in buffalo is still almost exclusively natural for the physiological problems in the application of the AI (Barile, Reference Barile2005; Drost, Reference Drost2007), we cannot exclude an inbreeding effect often due to the use of the same tested bulls.
Genotyping data, allele frequency, absolute (relative) frequencies of buffalo cows, lactations and tests across genotypes of the g.133 A > C SNP at the SCD promoter gene of the in Mediterranean river buffalo population
Association between SCD genotype and milk traits
The SCD genotype was significantly associated with milk yield (P = 0·02) (Table 2). In particular, buffalo cows with heterozygous genotype AC at the promoter of SCD locus showed the highest daily milk yield, with more than 2 kg/d compared with CC buffaloes. Such a difference accounts for about 28% more milk per day. On the contrary homozygous AA were slightly lower than AC. The behaviour of the three genotypes tended to remain constant throughout the whole lactation, as can be seen from the estimated lactation curves of the three genotypes (Fig. 3). A similar pattern has been observed also in other association studies involving polymorphisms at the SCD gene in Italian Holstein (Macciotta et al. Reference Macciotta, Mele, Conte, Serra, Cassandro, Dal Zotto, Cappio, Pagnacco and Secchiari2008) and at the OXT gene in river buffalo (Pauciullo et al. Reference Pauciullo, Cosenza, Steri, Coletta, Jemma, Feligini, Di Berardino, Macciotta and Ramunno2011).
Lactation curves of different SCD genotype for milk yield (kg/d)
Least squares means of milk yield (kg/d) for the three genotypes at the g.133A > C SNP in the promoter of SCD gene estimated with model (1)
a,bMeans within columns without a common superscript differ (Bonferroni adjusted P<0·05)
The allele substitution effect of the adenine into cytosine was about −1 kg/d (P < 0·01), and the contribution of the SCD polymorphism to the total phenotypic variance was 12% (Table 3). Such an effect is larger than the one reported for the SCD (Macciotta et al. Reference Macciotta, Mele, Conte, Serra, Cassandro, Dal Zotto, Cappio, Pagnacco and Secchiari2008) and also higher than the contribution reported for DGAT1 (Grisart et al. Reference Grisart, Coppieters, Farnir, Karim, Ford, Berzi, Cambisano, Mni, Reid, Simon, Spelman, Georges and Snell2002) on milk yield in dairy cattle.
Mean±se of milk yield (kg/d) for the substitution effect g.133A > C in the promoter region of the river buffalo SCD gene and contribution of this polymorphism to the phenotypic variance
† α: Substitution effect
‡ d: dominance effect
§ σ2: variance components associated to the genotype (SCD); to the individual buffalo cow (c), to residuals (e)
¶ r 2: contributions of genotype (SCD) and of individual buffalo cow (c) to the total phenotypic variance
The large effect of dominance on milk yield observed in the present work, more than 1·2 kg (P < 0·02), is also very interesting (Table 3). This offers a possible explanation of the over-dominance effect of the heterozygous (AC) on the best homozygous phenotype (AA). Often such an effect is not detected or considered to be not relevant because numerically much lower than the additive effect. However, it might have an impact on allele substitution effect in the population as recently reported in dairy cattle (Kuehn et al. Reference Kuehn, Edel, Weikard and Thaller2007). Moreover, lactation curves of different genotypes showed that the superiority of the heterozygous over the AA genotype is concentrated mainly around the lactation peak and, to a lesser extent, between 120 and 150 DIM (Fig. 3).
In previous research carried out on the same group of animals, Pauciullo et al. (Reference Pauciullo, Cosenza, D'Avino, Colimoro, Nicodemo, Coletta, Feligini, Marchitelli, Di Berardino and Ramunno2010) showed that CC genotype is associated to higher desaturation index in milk fatty acids, probably as consequence of an additional SP1 transcription factor in the promoter. It is also well known that the lipid metabolism and the de-novo biosynthesis of fatty acids are complex pathways and they are energetically very expensive. As suggested for dairy cattle (Kay et al. Reference Kay, Weber, Moore, Bauman, Hansen, Chester-Jones, Crooker and Baumgard2005; Macciotta et al. Reference Macciotta, Mele, Conte, Serra, Cassandro, Dal Zotto, Cappio, Pagnacco and Secchiari2008), in cows with greater desaturase activity, fewer nutrients are directed toward milk yield. In river buffalo, the CC genotype at SCD locus showed an increased Δ9-desaturase activity and higher milk monounsaturated FA content (Pauciullo et al. Reference Pauciullo, Cosenza, D'Avino, Colimoro, Nicodemo, Coletta, Feligini, Marchitelli, Di Berardino and Ramunno2010), but it is also characterized by lower milk yield. Therefore, our findings would seem to confirm this relationship.
Environmental factors affecting milk traits
Parity and year × calving season significantly affected buffalo daily milk yield (Table 4). Milk production tended to increase from first to later parities, reaching the maximum at the fourth calving. These findings are in agreement with previous reports on buffaloes (Catillo et al. Reference Catillo, Macciotta, Carretta and Cappio-Borlino2002; Dang et al. Reference Dang, Mukherjee, Kapila, Mohanty, Kapila and Prasad2010). Across all years, highest daily milk yields were observed for buffaloes calving in autumn whereas lowest values were for summer calvings (not reported for brevity). The seasonal effects on productive performances of buffaloes are not news. As reported in Mediterranean and Pakistan Nili–Ravi breed (Catillo et al. Reference Catillo, Macciotta, Carretta and Cappio-Borlino2002; Afzal et al. Reference Afzal, Anwar and Mirza2007), a depressive effect of high temperatures at the beginning of the lactation is recorded for calving in summer.
Least squares means of milk yield (kg/d) for the different levels of parity estimated with model (1)
A,B,C,DMeans within columns with different superscripts differ (P < 0·01)
Conclusions
Stearoyl-CoA desaturase genes have been intensively investigated in ruminants mainly through association studies between SNP and the fatty acid spectrum. This study reports an association between the g.133A > C SNP previously found in the promoter region of SCD gene and daily milk yield in the Italian river buffalo.
Little information is available for this species about association studies but also on the effectiveness of traditional selection schemes. On the other hand, much work has been done to improve recording, health, feeding and livestock systems.
This study represents one of the first indications of an association between a trait of economic importance and a candidate locus in river buffalo. A full characterization of the SCD promoter region was proposed, genotyping information for the SNP g.133A > C was provided and a quick method for allelic discrimination was set up. A significant association with daily milk yield has been found at SCD gene. The genotype AC showed an over-dominance effect with an average daily milk yield approximately 1·2 kg/d higher. The effect of the genotype was constant across lactation stages. Although such results need to be confirmed with large-scale studies in the same and other buffalo populations, they might be of great economic interest for the buffalo dairy industry.
This work was financially supported by the Italian Ministry for Agriculture and Forestry Policy—MiPAAF (INNOVAGEN project).
