Polymorphism of the STAT5A, MTNR1A and TNFα genes and their effect on dairy production in Bubalus bubalis

Abstract

The water buffalo is a fundamental resource, especially in developing countries, however, differently from other species, its genetic potential is still poorly investigated. In this work, we performed a candidate gene association study for milk composition in 491 female buffaloes. Animals were from four farms located in Southern Italy, where the Out-of-Breeding-Season-Mating technique is usually performed. We analysed three genes: (1) the signal transducer and activator of transcription 5A (STAT5A), (2) the tumour necrosis factor alpha (TNFα) and (3) the melatonin receptor 1A (MTNR1A). We confirmed the mutation at the MTNR1A gene and we found five novel single nucleotide polymorphisms (SNPs): one in the TNFα and four in the STAT5A. No associations were found for the SNPs in the MTNR1A and TNFα genes, while we identified a marked association with milk protein % for a C > T substitution at the STAT5A gene. At this locus, the TT buffaloes showed significantly higher protein percentage in milk. Conversely, this genotype class was the less frequent in the population. Moreover, an A > G substitution at the STAT5A showed an influence on reproductive seasonality, with the advantageous allele most frequent in the population, suggesting a possible effect of selection for this trait. The C > T substitution on STAT5A detected in present study could be used in marker assisted selection of Mediterranean Italian buffalo, and should be monitored to understand the reasons behind the low frequency of the favourable genotype at this locus and to stop this unfavourable trend in the population.

Keywords:

• River buffalo
• STAT5A
• TNFα
• MTNR1A
• single nucleotide polymorphism

Introduction

Water buffaloes (Bubalus bubalis) provide 13% of the total world milk production and, in Europe, 93% of animals are in Italy (faostat.fao.org). Here, the Mediterranean Italian buffalo breed is reared for milk production, mostly processed into mozzarella cheese. Advances in management techniques (Barile 2005) together with genetic selection (Iamartino et al. 2013), led to the high productive standards of these animals, which are exported to many countries. However, genetic improvement is still poorly exploited in this species, also due to the lack of phenotypes. Late age at maturity, long calving intervals and reproductive seasonality are other obstacles to the genetic improvement in the species (Zicarelli 1997, 2004; Campanile et al. 2005). During the past few decades, advances in molecular genetics have led to the identification of genetic markers associated with genes that affect traits of interest in livestock (Dekkers 2004). Molecular markers increased the accuracy of selection, especially at young ages, leading to a reduction in optimum generation length and an increase in rate of genetic gain (Goddard & Hayes 2002). Genomic prediction of breeding values is increasingly widely used in breeding programmes of many livestock species (Hayes et al. 2014, WCGALP), but buffalo genomic data are still lacking. In a recent work, authors demonstrated that the cattle SNP chip does not offer an optimal coverage of buffalo genome (Camargo et al. 2015). Hence, the candidate gene association study is still a powerful method in this species. Moreover, to identify the causative mutations within genes or regulatory sequences, SNPs within genes provide a higher power for association analysis, compared with using a SNP set with uniform genome-wide distribution (Jorgenson & Witte 2006).

In this work, we investigated the effects of several important genes on dairy and reproductive traits by genotyping a representative sample of Mediterranean Italian Buffalo, aiming to find polymorphisms that could be useful for the genetic improvement in this species. We analysed the following genes: signal transducer and activator of transcription 5A (STAT5A), tumour necrosis factor alpha (TNFα) and the melatonin receptor 1A (MTNR1A) subtype.

The STAT5A is a member of the interferon-τ (IFN-τ) and placental lactogen (PL) signal transduction pathway, which is very important in both milk production and fertility. The STAT proteins are widely studied for their role in cytokine signalling pathways, as they are transcription factors that act as signal transducers in the cytoplasm and transcription activators in the nucleus (Kisseleva et al. 2002). One single nucleotide polymorphism (SNP) in the STAT5A gene has been found to be associated with milk protein and fat percentages and with embryonic survival (Khatib et al. 2008; He et al. 2012). An effect of this gene on milk yield was also found in goat (An et al. 2013). The TNFα gene affects both gonadotrophin production and male reproductive function in mouse, man and dairy cow (Tronchon et al. 2008; Zalata et al. 2013). The MTNR1A, that affects seasonal reproductive activity in sheep, goat and buffalo (Pelletier et al. 2000; Carcangiu et al. 2009, 2011), has been associated with milk protein percentage in buffalo (Zetouni et al. 2014).

Materials and methods

Animals

SNP detection was conducted on 12 animals by a sequence alignment. SNP genotyping and association analysis were performed on 491 female river buffaloes sampled from four farms in Southern Italy, with a mean of 114 animals per farm.

Blood samples analysed in present work were collected by a veterinarian within the routine tests carried out at the farms. Sampling was performed with the consent of the breeders and in compliance with relevant national and international regulations on animal welfare requirements (Eurpean Council 1986).

DNA extraction

DNA extraction was carried out from 200 μL of frozen whole blood using a commercial kit (Promega ReliaPrep Blood gDNA Miniprep System, Madison, WI), according to the manufacturer’s instructions. The concentrations measured through agarose gel were around 50 ng/μL.

SNP detection and genotyping

SNP detection was performed by PCR amplification and sequencing. Amplified regions, primers and annealing temperatures are reported in Table 1. The PCR reaction mix (20 μL) comprised: 1 μL of gDNA, 5× PCR Buffer (Promega, Madison, WI), 5 mM MgCl₂, 0.4 μL of each primer, dNTPs each at 5 mM, 0.2 μL of Taq DNA Polymerase (Promega, Madison, WI). PCR products were purified and sequenced in outsourcing (BMR Genomics; http://www.bmr-genomics.it/). The alignment of sequences was carried out using the software BioEdit (Carlsbad, CA) (Hall 1999).

Table 1. Amplified regions, primers and annealing temperatures used for detecting single nucleotide polymorphisms in present work.

Gene	Amplified region	Sense primer 5′	Antisense primer 5′	PCR product	Annealingtemperature
MTNR1A	Exon 2	GTGGTGAGCCTGGCAGTT	ATGGAGAGGGTTTGCGTTTA	611 bp	65 °C TD^a
TNFα	Exons 3–4	CAAGTAACAAGCCGGTAGCC	GGACACCTTGACCTCCTGAA	741 bp	65 °C TD
STAT5A	Exons 2–3	GGACATGGCGGGCTGGAT	GGAGCTGCGTGGCATAGT	318 bp	60 °C TD
	Exons 8–9	GTGTGAGAAGTTGGCGGAGA	TTGCGGGTGTTCTCGTTCT	740 bp	65 °C TD
	Exons 10–11	TGAGTGCAGCGGGGAGAT	CTTCACCTGGAACACAAGCTC	305 bp	65 °C TD
	Exons 12–13	GTCCCTTCCCGTGGTTGT	GGTTGAACTGGGACCAGGA	156 bp	65 °C TD
	Exons 14–15	ACTTGCCCGGCTGGAACTAC	AGAGTCAAACTTCCAGGCGATG	108 bp	65 °C TD
	Exon 16	CTGACCGTAACCTGTGGAATC	CAAGCACAGGAGTGTAGTACTTGG	285 bp	65 °C TD
	Exon 18	GTTTGTGAGCGCCTCTGC	TTCTGTGGGTACATGTTATAGTGAGG	1296 bp	64 °C
	Exon 19	CCTGACCCGGTGCTCGAC	GTGTACATGGGCTGCCTGCAA	817 bp	62 °C

a Touchdown (TD), in all TD PCR performed, starting temperature was 5 °C higher in respect of the optimum ones, 10 cycles with decreasing temperature of 0.5 °C each.

Genotyping of the SNPs was performed on genomic DNA extracted from whole blood samples. The analysis was carried out in outsourcing, using the KASPar SNP genotyping system developed and patented by LGC Genomics (www.lgcgenomics.com).

Phenotypes

The Italian Buffalo Breeders’ Association (ANASB) provided the official data of the milking records of the animals genotyped. For the association analysis, the following parameters were used: farm, birth date, date of calving, lactation number and full-lactation (270-day) milk production traits (milk yield kg, protein kg and % and fat kg and %). These phenotypes were provided for 391 out of 491 buffaloes genotyped, as for 100 buffaloes only the date of the first calving was available. After checking quality of the phenotypic data, we used a total of 1064 records on 391 buffaloes to evaluate production traits with respect to the genotypes. For fertility evaluation, calving intervals were calculated as the distance in days between two consecutive dates of calving. The same procedure was applied to every animal having more than one calving date recorded. Since under optimal conditions buffalo can resume ovarian activity after calving by 30–90 days (Moioli et al. 1998), we considered only calving intervals falling in a range between 330 and 1000 days, for a total of 946 observations on 358 animals. The dates of calving were also used to consider the reproductive seasonality. The month of calving was present for all the recorded lactations and for all buffaloes, for a total of 1441 observations on 491 animals. In the analysed farms, breeders apply the Out-of-Breeding-Season-Mating (OBSM) technique, which consists in suspending sexual promiscuity in the herd during autumn to allow milk production in spring. Subsequently, the SNP association analysis was performed between ‘seasonal’ (from August to December) and ‘out of season’ (from January to July) dates of calving.

Statistical analysis

Population genetic parameters were calculated using PowerMarker v3.25 (Liu & Muse 2005). Minor allele frequency, expected and observed heterozygosity, and Hardy–Weinberg (HW) equilibrium exact p values were calculated. Associations between the phenotypes and the polymorphism analysed were tested using JMP^®, Version 12 Pro (SAS Institute Inc., Cary, NC, 1989–2007) with the following mixed linear model: (1) where Y_ijklm = calving interval (days), seasonality (seasonal, out of season), milk yield (kg), protein or fat yield (kg) and protein and fat content (%); μ = overall mean; FARM_i = fixed effect of the farm (four levels); bCI  =  covariable represented by the calving interval; NLACT_j = fixed effect of the number of lactation (from 1 to 7); SEA_k = fixed effect of the calving season (two levels), excluded from the model when used as a variable; B_l = random effect of the animals (N. 391); SNP_m = fixed effect of the SNP genotype (three levels); and e =  random residual. It was not possible to use the effect of the sire because this information was not provided. This is a critical point in buffalo breeding. Rates of sire and dam misidentification amount to 24% and 20% in the Italian buffalo population (Parlato & Van Vleck 2012). Average gene substitution effect (α) was estimated using the same model, but with the gene effect treated as a covariate (Banos et al. 2008; Pauciullo et al. 2012). The coding of the three genotypes (0, 1, 2) was based on the number of copies of the first allele in alphabetical order (Dagnachew et al. 2011). The additive effect (a) was estimated by comparison of the means of the traits value for homozygote, i.e. a = 1/2(BB − AA) and the dominance effect (d) as AB −1/2(AA + BB). The ratio |d/a| was considered to indicate gene effects (Stuber et al. 1987): |d/a| < 0.2, additive; 0.2 < |d/a| < 0.8, partial dominance; 0.8 < |d/a |< 1.2, dominance; |d/a| > 1.2, overdominance (Fontanesi et al. 2012). In order to isolate the interaction of fixed effects, we evaluate the interactions between the fixed factors considered (farm-seasonality, farm-parity, parity-seasonality).

Results

We identified six novel SNPs, by aligning the PCR fragments obtained (Table 2). The SNP STAT5A_c128plus179G was not in HW equilibrium, with significantly different values for observed and expected heterozygosity, and the MTNR1Ac.318C > T was close to the threshold of significance (Table 3). No associations were found for the polymorphisms in the MTNR1A and the TNFα genes. The SNPs found at STAT5A gene are not in linkage disequilibrium.

Table 2. Characteristics of six SNPs found in Bubalus bubalis MTNR1A, TNFα and STAT5A genes, in a sample of 12 animals collected in South of Italy.

Gene	SNP name	Location^a	SNP type	aa. change	References	RefSeq-SNPIdentification^b
MTNR1A	c.318C >T	Exon 2	C/T	No	Carcangiu et al. (2011)	GU817415
TNFα	c.323G > A	Exon 4	A/C	No	Present work	2019323453
STAT5A	c.128 + 179C	intr2–3	G/C	–	Present work	2019323457
	c.924A > G	ex8	C/T	no	Present work	2019323456
	c.989 + 344A	int8–9	C/T	–	Present work	2019323455
	c.1342 + 99T	intr10–11	A/G	–	Present work	2019323454

a Based on the bovine sequence (gene ID: 282375).

b Reference sequence (RefSeq) and SNP identification (ss number) in the National Center for Biotechnology Information (NCBI) database (Pruitt et al. 2014).

Table 3. Population genetics parameters: genotype frequencies, minor allele frequencies (MAF), observed and expected heterozygosity (Het obs and Het exp, respectively), and Hardy–Weinberg (HW) equilibrium exact p values at five SNPs found in Bubalus bubalis MTNR1A, TNFα and STAT5A genes, genotyped on 491 total Mediterranean Italian buffaloes.

SNP	No. of obs.	Genotype and frequency			Allele and MAF	Het Exp	Het Obs	HW exact p Value
MTNR1A_c318CT		CC	CT	TT	T
	487	0.520	0.378	0.103	0.292	0.413	0.378	.050
TNFα_c323CA		AA	AC	CC	A
	491	0.0163	0.1894	0.7943	0.111	0.197	0.189	.364
STAT5A_c128plus179G		CC	CG	GG	C
	467	0.0728	0.4925	0.4347	0.319	0.435	0.493	.004
STAT5A_c924CT		CC	CT	TT	T
	483	0.4224	0.4389	0.1387	0.358	0.460	0.439	.273
STAT5A_c989plus344C		CC	CT	TT	T
	485	0.4742	0.4103	0.1155	0.321	0.436	0.410	.172
STAT5A_c1342plus99A		AA	AG	GG	G
	484	0.5496	0.3678	0.0826	0.267	0.391	0.368	.199

The STAT5A_c989plus344C locus was associated with milk protein percentage (p < .01) showing significantly higher performances for the TT genotype (Table 4). The fixed factors of farm and number of lactation were statistically significant in the model, with p < .01 and p < .05, respectively. The interactions of fixed effects were not statistically significant. The mean value for the protein percentage parameter was 4.61 in the analysed population. Even if the allelic substitution effect, equal to 0.007, was not statistically significant, the C allele showed a negative overdominance effect (d/a = –5). The genetic additive effect was equal to 0.01 and the SNP explains 2.4% of the total variability of our data with the model used. The TT buffaloes at this locus showed also a significantly reduced reproductive seasonality (p < .001), as shown in Table 5. The SNP at position c.1,342 + 99A affected reproductive seasonality as well (p < .05), even if the latter association was not significant considering the Bonferroni-adjusted significance level of 0.0125.

Table 4. Least squares means and standard errors (S.E.) of the genotypes at the SNP STAT5A_c989plus344C locus significantly associated with protein percentage (PP) in a Mediterranean Italian Buffalo population, estimated by a mixed linear model.

N. obs^a	N. animals	LSM ± S.E.	N. obs.	N. animals	LSM ± S.E.	N. obs.	N. animals	LSM ± S.E.
C:C			T:C			T:T
279	121	4.60AB ±0.01	273	111	4.56B ± 0.01	76	28	4.62A ± 0.02

A, B means with different superscripts differ (p < .01).

a Number of observations.

Table 5. Percentage of ‘seasonal’ (from August to December) and ‘out of season’ (from January to July) calvings for the different genotypes at the STAT5A_c.989 + 344C and STAT5A_c.1342 + 99A loci in a Mediterranean Italian Buffalo population, estimated by a mixed linear model.

Locus	p	Genotype	N. animals.	N. obs.	Out of season, %	Seasonal, %
STAT5A_c989plus344C	.0006	C:C	159	395	30.2	14.6
		T:C	145	384	28.0	15.5
		T:T	34	103	9.2	2.5
STAT5A_c1342plus99A	.0403	A:A	510	183	40.4	16.8
		G:A	327	128	22.9	13.8
		G:G	53	27	3.7	2.2

We did not confirm the influence of the MTNR1Ac.318C > T polymorphism on reproductive seasonality.

Discussion

In this work, we detected five novel SNPs on candidate genes for dairy traits in river buffalo. The SNP STAT5A_c128plus179G showed an observed heterozygosity higher than the expected, with the homozygous class for the minor allele C highly under-represented in the population. This can be due to a selection against this allele, which we found in the population mostly in heterozygous form. Other two loci at the STAT5A were associated with the parameters analysed: the STAT5Ac.989 + 344C located between the 8th and 9th exon, while the STAT5Ac.1,342 + 99A between the 10th and the 11th. It is well established that non-coding SNPs may affect important mechanisms such as transcription, translation and splicing. It is now well established the potential effect of pre-mRNA splicing on phenotype traits (Nissim-Rafinia & Kerem 2002; Ho et al. 2011). Moreover, non-coding variations that are close to a functional mutation are almost always inherited together with the functional mutation itself (Kwon & Goate 2000).

In this work for the first time, we assessed the importance of the STAT5A gene on dairy parameters in buffalo (Bubalus bubalis) through evidence of a significant association of the SNP STAT5Ac.989 + 344C with the protein percentage in milk. This is a fundamental trait for mozzarella cheese production, which is the main income from buffalo in Italy, the cheese being exported worldwide. The effect is evident by the significantly higher performances of the animals carrying the TT genotype. The major allele C showed an overdominant effect, with a milk protein content of the TC genotype lower than the two homozygous, confirmed by a |d/a| ratio higher than the threshold value of 1.2 (Stuber et al. 1987). This may be the reason why the allelic substitution effect was not statistically significant, as the overdominance effect could limit the ability of the infinitesimal model to estimate the additive genetic effect of the SNP.

The STAT5A protein was identified in the mammary gland as a regulator of milk protein gene expression (Watson 2001) and is a member of the IFN-τ and PL signalling pathway. It is involved in signal transduction within many cells, including the uterus and mammary epithelial cells. Our results in buffalo confirm what was observed for the STAT5A gene in other species, such as cattle (Khatib et al. 2008; Bao et al. 2010; Oikonomou et al. 2011; He et al. 2012) and goat (An et al. 2013). In addition, another member of the STAT family, the STAT1 gene, has been recently found to be associated with milk production traits in buffalo (Deng et al. 2016).

The STAT5Ac.989 + 344C, together with the c.1,342plus99A locus, showed significant effects on reproductive seasonality, even though the reproduction technique used in the analysed herds could have influenced the effect. Given the association of this genotype with higher protein percentage in milk, further analyses on a population where the OBSM technique is not applied would be interesting to confirm the trend highlighted in present work. Moreover, the most effective genotypes for milk production seem to have a lower sensitiveness to the photoperiod. The same trend is observed for the c.1,342plus99A polymorphism, where the most represented genotypes showed more than 74% of the calvings in the out of breeding season period, even if considering the most conservative Bonferroni-adjusted significance level. This may indicate a selection of the less seasonal animals. In fact, all the analysed farms apply the OBMS technique to allow production in spring, the most favourable period for market demand of mozzarella cheese. The efficiency of this practice is influenced by the management level of the farm. As one of the four farms in the study showed significantly higher percentage of animals calving in the out of breeding season period, we repeated the analysis excluding all the observations of this farm. Both SNPs were still significant after Bonferroni’s correction (p < .0001). In all cases, the SNP on MTNR1A did not show any influence on reproductive seasonality. Once again it is worth mentioning that the conditioning of reproduction performed in the farms could have mitigated the effects of this SNP.

The role of STAT5A gene in reproduction is well known in ruminants. In particular, this gene influences embryonic development and the signal transduction pathway of IFN-τ, which has a key role in the initiation and maintenance of pregnancy (Spencer & Bazer 2004). Further analyses, together with the collection of more phenotypes for fertility evaluation, could assess the trend highlighted in the present work on seasonality and unveil other interesting effects of the STAT5A polymorphisms in buffalo also on the reproduction and seasonality. Moreover, the availability of reliable phenotypes on oestrus behaviour could evidence, besides the effect on production, further interesting associations with oestrous expression, as already found in cattle for the STAT5A gene (Homer et al. 2013). This would be remarkable in buffalo species, where the difficulty in oestrus detection is a limiting factor for the use of artificial insemination.

Conclusions

In this work, we demonstrated the importance of the STAT5A gene in buffalo species, by assessing an association with protein content in milk. This trait is fundamental for mozzarella cheese production, which is the main product of buffalo breeding in Italy exported in many countries. The C > T substitution at the STAT5A found in present study could be used in marker assisted selection of Mediterranean Italian buffalo, by increasing its frequency in the population. However, an accurate recording of phenotypic parameters, as well as genealogies registration, are critical elements that must be improved to increase the genetic progress in buffalo. We also highlighted a possible influence on the reproductive seasonality, which deserves to be better evaluated in a population where the conditioning of reproduction, usually applied in Italian farms, is not performed. No associations were found for MTNR1A and TNFα genes with the analysed parameters. The availability of a high-throughput SNP platform in the next future for buffalo species will certainly enhance the perspectives for genetic improvement of this species.

Disclosure statement

No potential conflict of interest was reported by the authors.

Additional information

Funding

This research was supported by the CISIA Project – Integrated Knowledge for Sustainability and Innovation of Made in Italy Agribusiness – VARIGEAV, CNR-ISPAAM (National Research Council – Institute for Animal Production Systems in Mediterranean Environments).