Association analysis of polymorphism in thyroglobulin gene promoter with milk production traits in riverine buffalo (Bubalus bubalis).

Polymorphism within the promoter region of bovine thyroglobulin has been reported to be associated with milk and meat quality. In this study, we investigated the genetic variation within thyroglobulin promoter region of swamp and riverine buffaloes using PCR-SSCP technique and sequencing, and also analyzing association of polymorphism with the milk production traits. The study revealed four conformational patterns, A, B, C, and D among 323 buffaloes of two riverine breeds and different swamp populations. The frequency of SSCP variant 'A' was found to be invariably high among all buffalo populations. Variant 'C' was found to be absent in pure swamp population and present with higher frequency among riverine dairy breeds Mehsana and Nili Ravi. Frequency of D variant was observed to be highest in buffalo population, representing riverine and hybrid types. Sequencing of three representative PCR products of each of the SSCP variants, revealed three polymorphic sites responsible, 33C > T, 176G > T and 221C > T, in the buffalo TG promoter region. Further, association studies of SSCP variants with various milk production and milk quality traits indicated significant effect on fat percentage in buffaloes belonging to Mehsana and Nili Ravi dairy breeds. The preliminary results also showed the substantial variations in the distribution of SSCP variants' frequencies across swamp and riverine buffaloes, two distinct populations being reared for meat and milk production, respectively.

Introduction

Buffalo is an important livestock species contributing to meat, milk and draught power in Southeast Asia. Fat percentage is important parameter to assess the quality of milk as well as meat of both cattle and buffalo. A QTL region associated with fat deposition in skeletal muscles, has been identified in the centromeric region of bovine chromosome 14 (BTA14), among multiple cattle populations (Moore et al., 2003, Casas et al., 2005). Genes lying in this region include diacylglycerol-O-acyltransferase (DGAT1), thyroglobulin (TG) and adipose fatty acid binding protein (FABP4), found to be associated with fat percentage in both beef and dairy cattle (Barendse, 1999, Michal et al., 2006). Thyroid hormones play an important role in regulating the metabolism and can affect adipocyte growth, differentiation, and homeostasis of fat depots. Among these genes, TG is a glycoprotein hormone, synthesized in thyroid follicular cells and carrier for both triiodothyronine (T3) and thyroxine (T4), stored in the thyroid gland. Genetic variation in TG has been associated with back fat thickness and marbling in beef cattle as well as milk traits in dairy cattle (Moore et al., 2003, Barendse et al., 2004, Rincker et al., 2006, VanEenennaam et al., 2007). The genetic variations responsible are located in the 5′promoter region and 3′UTR of the TG gene, which have been widely used in marker assisted selection (MAS) programs to improve the predictability of marbling level and eating quality in beef (Barendse, 1999). An allele of the TG gene has also been identified to be having a significant association with marbling score in cattle (Gan et al., 2008). In Japanese Black cattle, marbling has been associated with both T3 and T4 (Mears et al., 2001). The polymorphism C/T (PsuI-RFLP) in 5′ upstream region was associated with higher marbling scores in cattle that had the homozygous cytosine (C) allele as favourable one (Wood et al., 2006).

India has two different populations of water buffaloes — swamp and riverine types, differing in their chromosome numbers as well as in habitat and utility. Phenotypically also two buffalo populations differ, with swamp buffaloes mainly found in Northeast region of India, being short and stout are utilized for draught and meat purpose, whereas riverine buffaloes not so muscular in body structure, are reared primarily for milk. Earlier reports have also shown significant variations in the allelic frequencies of polymorphic loci in genes governing immune response and other traits across swamp and riverine buffaloes (Dubey et al., 2013), indicating two populations with variable genetic structure.

Majority of Indian buffalo population is riverine type, contributing more than 50% to total milk production. Buffalo species is well known for high fat percentage in its milk, quality of meat low in cholesterol and candidate gene polymorphism has been exploited to find the association with these production traits (Tanpure et al., 2012). Due to the important role of thyroglobulin in fat metabolism and association of polymorphism with meat and milk quality reported in cattle and because of the fact that no such analysis has been carried out in closely reported buffalo species the present study was undertaken to identify polymorphism within promoter region of thyroglobulin gene across riverine and swamp buffaloes and analyze its association with milk production and fat percentage traits in buffalo.

Materials and methods

Animals, blood sample collection and DNA isolation

Blood samples were collected from unrelated buffaloes belonging to different breeds/populations from their native breeding tracts of different agro-climatic regions of the country. For SSCP analysis, the breeds/populations included were riverine type dairy breeds — Mehsana (136) and Nili Ravi (82) along with swamp type of Northeast Indian states including, Assam-riverine or hybrid (48), Assam – pure swamp (41) and Mizoram – pure swamp (16), bordering China and Myanmar (Table 1). Since buffaloes in India are mostly dairy type, riverine breeds of Mehsana (N = 144) and Nili Ravi (N = 22), for which milk production and fat records existed, were included for association analysis studies. Genomic DNA was extracted from the whole blood using standard SDS-Proteinase-K digestion and phenol/chloroform extraction procedure (Sambrook and Russel, 2001). Quality and quantity of DNA was assessed by measuring OD260 and OD280 using a Nanodrop, ND-1000 UV–vis Spectrophotometer (Thermo Scientific) as well as by running on ethidium bromide stained 0.8% agarose gel. The study was conducted following the ethical and regulatory guidelines in place at the institute.

Table 1

Frequency of SSCP variants in thyroglobulin promoter in different riverine and swamp buffalo breeds/populations.

Breed/population	Total no.	Variants (no. of observations)				Frequency
A	B	C	D	A	B	C	D
Mehsana	136	94	19	17	6	0.69	0.14	0.13	0.04
Nili Ravi	82	68	3	10	1	0.83	0.04	0.12	0.01
Total	218	162	22	27	7	0.75	0.10	0.12	0.03
Assamese (riverine/hybrid)	48	34	5	3	6	0.71	0.10	0.06	0.13
Swamp	57	39	5	0	13	0.69	0.08	0.00	0.23
Total	105	73	10	3	19	0.69	0.10	0.03	0.18
Overall	323	235	32	30	26	0.72	0.10	0.10	0.08

Single strand conformation polymorphism analysis and sequencing

A set of oligonucleotide primers (forward: 5′-AGAGGGGAAAGGGGGATGACT-3′ and reverse: 5′-GGGGGTGTGGCTTCATAATTC-3′) was designed from the cattle thyroglobulin gene promoter sequence using PrimerSelect program of Lasergene software (DNASTAR Inc., Madison, WI, USA), to amplified 300 bp region, comprising part of promoter region, − 423 bp upstream to exon 1 start site (Dubey et al., 2014). Polymerase chain reaction was performed in a final reaction volume of 25 μL, containing ~ 100 ng of genomic DNA, 5 pmole of each primer, 200 μM of each dNTPs, 2.5 μL of 10 × buffer with 15 mM MgCl2 and one unit of Taq DNA polymerase (Bangalore Genei, India). Amplification was performed with an initial denaturation at 95 °C for 2.5 min followed by 35 cycles of 94 °C for 30 sec, annealing temperature at 58 °C for 30 sec and extension at 72 °C for 1 min, with a final extension of 5 min at 72 °C.

The PCR products were visualized on 2% ethidium bromide stained agarose gel. Amplified PCR products were subjected to single-strand conformational polymorphism (SSCP) analysis after optimization of non-denaturing PAGE conditions. The electrophoresis was carried out in 8% PAGE gel by using Protean II vertical electrophoresis unit (Bio-Rad, USA) using 1 × TBE buffer. Gels were silver stained (Sambrook and Russel, 2001), dried and scored manually for SSCP variants. The SSCP variants were further sequenced for the identification of single nucleotide polymorphisms within promoter region of thyroglobulin gene using the PCR products from three representative samples of each variant, representing Mehsana, Nili Ravi riverine and swamp buffalo populations. Purified PCR products were sequenced using BigDye terminator cycle sequencing kit (Applied Biosystems, USA) on ABI 3100 Genetic Analyzer and the raw sequence data was edited manually using Chromas Ver. 1.45 (http://www.technelysium.com. au/chromas.html). Each of the identified SSCP variants was sequenced from both ends using both forward and reverse primers and submitted to NCBI, GenBank (Table 2). Multiple sequence alignment of the buffalo sequences along with cattle sequence was carried out using MegAlign program of Lasergene software (DNASTAR, Inc., Madison, WI, USA) and different transcription binding factors were analyzed and compared using GPMiner and p-Match software (http://gpminer.mbc.nctu.edu.tw/, http://www.gene-regulation.com/pub/programs.html), as described previously for buffalo (Dubey et al., 2014).

Table 2

Nucleotide changes observed at three loci in different SSCP variants of bubaline thyroglobulin promoter and their GenBank accession numbers.

SSCP variants	Nucleotide position			Accession numbers
33	176	221
Variant A	C	G	C	KC607819
Variant B	C/T	G	C	KC607820
Variant C	C	G/T	C/T	KC355364
Variant D	C/T	G/T	C/T	KC355365

Statistical analysis

For association analysis of SSCP variants with dairy performance traits, data related to different animals belonging to the Nili Ravi and Mehsana breeds of riverine buffalo, were collected and classified according to different seasons, years and herds etc. Effects of different non-genetic factors like seasons (rainy, mid-June to September; winter, October to February; summer, March to mid-June), year (2001–2005), parity and influence of different management practices in the farm and field conditions were included in the analysis. The effects of non-genetic factors were estimated by least squares analysis of variance for non-orthogonal data using a complete fixed model as described by Harvey (1987). The effect of different SSCP variants on various production traits milk yield, milk fat percentage and total fat yield at late stage of lactation, i.e. 305 days, was analyzed by least squares method using a mixed model with random effect of sire and fixed effect of SSCP variant or genotype as given below.

Y is the dependent trait under study (milk yield, fat percentage and fat yield at different stages of lactation of kth individual with jth genotype and ith sire), μ is the overall population mean, Sis the random effect of the ith sire, Gj is the fixed effect of the jth genotype, and e is the random error, assumed to be normally independently distributed with mean zero and constant variance, i.e., NID (0, σe2). Results were documented for significance at p-value 0.05.

Results and discussion

Distribution of SSCP variants across riverine and swamp buffaloes

SSCP analysis of 300 bp buffalo thyroglobulin promoter revealed presence of four variants, designated as A, B, C and D, in all the four breeds/populations of buffalo (Fig. 1). The overall frequency of the variant A was observed to be highest (0.72), followed by variant B (0.10), variant C (0.10) and variant D (0.08) respectively. The variant C was found to be absent in all the pure swamp populations of Northeast region under this study, whereas Assamese hybrid and riverine buffaloes of Northeast region had a frequency of 0.06 for variant C. The observed frequency of variant A was, 0.75 among different riverine and slightly lower, 0.69 in swamp breeds/populations. Variant D was observed more frequent among all the swamp, riverine and hybrid buffalo populations of Northeast region, whereas the frequency of variant B was observed similar (0.10) among riverine and swamp buffaloes (Table 1). Since allelic frequency analysis for the SSCP variants across riverine and swamp type buffaloes in this study showed higher frequency of variant D in swamp and variant C in riverine buffaloes comparatively, these variations could be responsible for the differences in the production traits, since swamp buffaloes are mainly reared for meat and draught, whereas riverine are used as dairy animals in India. This is an important finding, since in India riverine buffalo is a major dairy animal, contributing more than 50% of the total milk produced, mainly due to more liking of buffalo milk taste, having higher fat and protein percentage compared to cattle. Earlier reports have also indicated existence of substantial variation in the TG gene promoter region, across bovine species with some of the changes present within the important putative transcription factor binding sites and fixation of favourable allele for marbling at 422C > T locus in buffalo (Dubey et al., 2014). Taking lead from these earlier findings, the region further − 423 nucleotides upstream was targeted in the present study. Genetic variation within promoter and 3′UTR of TG gene has been widely used in marker assisted selection (MAS) programs to improve the marbling level and eating quality of beef cattle (Barendse, 1999).

Fig. 1

Non-denaturing polyacrylamide gels, showing different single-strand conformation polymorphism (SSCP) variants in the promoter region of buffalo thyroglobulin gene.

Sequence analysis of SSCP variants

Sequencing of representative SSCP variants, confirmed the nucleotide changes responsible for exhibiting different patterns in the gel. Sequences of amplified SSCP variants of buffalo thyroglobulin gene were compared with reported cattle sequence, revealing deletion of nucleotide A at position 52 (671 from exon 1, as per Accession number JX090181) in buffalo (Fig. 2). SSCP variants were found to be arising due to nucleotide polymorphisms, 33C > T, 176G > T and 221C > T (Table 2). Overall riverine buffalo sequence exhibited maximum homology with cattle and having similar transcription factors binding sites (Supplementary Table 1). However, deletion of single nucleotide A at the upstream position of 671 nucleotides from exon1 in buffalo (Dubey et al., 2014) as compared to cattle leads to generation of an important transcription factor binding site, SREBP1, reported to be playing an important role in lipid metabolism (Zhou et al., 2012). Further studies are also required to establish the impact of deletion within putative transcription factor SREBP1 binding site, on thyroglobulin gene expression in cattle and buffalo, two important dairy species of India, differing in the fat quality as well as fat contents in meat as well as milk. Two other SNPs identified in buffalo, 176G > T was present within another transcription factor binding site, c-Ets-2 and 221C > T was in the HSF transcription factor binding site in buffalo (Supplementary Table 1), could also affect the thyroglobulin expression and need to be further explored. All these SNPs identified are novel to buffalo, not reported previously.

Fig. 2

Multiple alignment of nucleotide sequences of different SSCP variants (variant A, variant B, variant C and variant D) of buffalo thyroglobulin gene promoter with reported cattle sequence of the same region. The position on top refers to cattle sequence, having one nucleotide insertion at position 52 compared to buffalo.

Association analysis

The SSCP variants identified were further analyzed in Mehsana and Nili Ravi buffaloes with production records, for association with milk production. Least square means (LSM) for milk traits linked to different SSCP variants obtained by statistical analysis revealed, mean total milk yield to be maximum for variant A (1897.61 kg) followed by variant C (1879.19 kg), B (1818.26 kg) and D (1615.92 kg) respectively (Supplementary Table 2). The 305 day milk yield (standard lactation length) also showed a similar trend without any significant variation. No significant association of variants with different traits was noticed, except for fat percentage, although differences among mean values of animals belonging to different SSCP variant groups were observed. There was significant (p < 0.05) effect of SSCP variants on both 305 days as well as total fat percentage (Table 3). Most of the other milk quality traits, like total fat yield, also showed non-significant association with SSCP variants. This could be due to large standard error observed because of small sample size for two (variant B-17 animals and variant D-7 animals) out of four variants. Total fat and 305 day fat percentages were found highest in variant B followed by variants D, A and C, respectively. In contrast, total fat yield was highest in variant B (130.38 kg), followed by variant A (126.70 kg), C (124.31 kg) and D (110.86 kg), respectively. Similar trend was observed with respect to 305 day fat yield in different SSCP variants.

Table 3

Analysis of variance (ANOVA) of different first lactation traits for herd, season and TG SSCP variant effects.

Trait	Source	DF	Mean squares	F-value	p-Value
Total milk yield	Herd	2	98,835.33	0.444	0.643NS
Season	2	116,391.40	0.522	0.594NS
SSCP variants	3	184,595.70	0.829	0.483NS
305 day milk yield	Herd	2	242,368.28	1.175	0.311NS
Season	2	114,312.89	0.554	0.576NS
SSCP variants	3	147,414.22	0.715	0.548NS
Lactation length	Herd	2	6834.92	3.729	0.026⁎
Season	2	2855.34	1.558	0.214NS
SSCP variants	3	3153.48	1.720	0.163NS
305 day fat percentage	Herd	1	0.04	0.092	0.762NS
Season	2	0.20	0.485	0.617NS
SSCP variants	3	1.66	3.959	0.010⁎
Total lactation fat percentage	Herd	1	0.04	0.102	0.750NS
Season	2	0.21	0.506	0.604NS
SSCP variants	3	1.66	3.986	0.009⁎
305 day fat yield	Herd	1	1067.51	0.888	0.348NS
Season	2	186.16	0.155	0.857NS
SSCP variants	3	681.25	0.567	0.642NS
Total fat yield	Herd	1	15.07	0.012	0.914NS
Season	2	176.31	0.138	0.871NS
SSCP variants	3	897.63	0.703	0.555NS

Values significant at p < 0.05. NS (Non-significant).

Previous studies in Jersey and Hungarian Simmental breeds of cattle have shown significant role of TG marker 422 C > T, in fat percentage and in 305 day milk production (Anton et al., 2012). Also similar results have been reported with CC animals showing significantly higher 305 day milk yield (kg) than TT cows. As for 305 day milk yield (kg), 305 day fat yield (kg), and maximum daily milk yield (kg) CC animals produced significantly higher values than other genotypes (Farkas et al., 2009). But Holstein sire family in another study did not show significant association between TG variants and milk production traits (Khatib et al., 2007). Previous studies in buffalo, whereas have shown fixation of favourable C allele at 422 locus (Dubey et al., 2014), indicating the need for novel efforts to assess the variability among the genes regulating important production traits in buffalo.

Preliminary association studies of genetic variation in the analyzed promoter region of buffalo have shown significant association with fat percentage of the milk quality trait in two dairy breeds of riverine buffalo. Since sample size is low, further investigations testing the variants' effects on large numbers of samples could help in finding the novel genetic variation to be responsible for meat as well as other milk production traits in buffalo, relatively an unexplored species.