Growth hormone (GH) is encoded by the GH gene belonging to the growth hormone-prolactin gene family: prolactin (PRL) and GH are hormones/cytokines, which coordinate a wide range of biological processes in vertebrates (Vijayakumar et al. Reference Vijayakumar, Novosyadlyy, Wu, Yakar and LeRoith2010). GH promotes somatic growth, influencing muscle and bone metabolism and lactation (Giustina & Veldhuis, Reference Giustina and Veldhuis1998) and, in sheep, it is involved in the response to stress (Carcangiu et al. Reference Carcangiu, Vacca, Parmeggiani, Mura, Pazzola, Dettori and Bini2008). The PRL and GH gene family is characterised by species-specific gene expansion probably due to gene duplication and natural selection (Otha, Reference Otha1993). The human genome contains a cluster of five GH genes, while cattle have only one (bovine GH, bGH), and sheep (ovine GH, oGH) (as well as goats) display a copy number variant (CNV) polymorphism, with the Gh1 allele carrying the only GH1 copy and the Gh2 duplicate allele carrying the GH2-N and GH2-Z gene copies, separated by a 3·5 kb long DNA intercopy region (Valinsky et al. Reference Valinsky, Shani and Gootwine1990). The sheep GH locus was mapped to chromosome 11q25 and has 5 exons spanning approximately 2 kb for the Gh1 allele and 7·5 kb for the duplicate Gh2 allele (Ofir & Gootwine, Reference Ofir and Gootwine1997). The products of the two sheep GH gene copies show different biological properties and differ in their receptor-binding ability (Reicher et al. Reference Reicher, Niv-Spector, Gertler and Gootwine2008) and in stage and organ of expression: GH1/GH2-N is expressed in the pituitary, while GH2-Z is expressed in the placenta (Lacroix et al. Reference Lacroix, Devinoy, Servely, Puissant and Kann1996; Gootwine, Reference Gootwine2004). Sheep GH1, GH2-N and GH2-Z genes share the same nucleotide sequence, except for nucleotide (nt) variations which do not represent a distinctive feature. The only element that allows distinguishing GH1/GH2-N from GH2-Z is a segment of the promoter region of GH2-Z, which is included in the intercopy region (Marques et al. Reference Marques, Santos, Carolino, Belo, Renaville and Cravador2006).
Genetic variants at the GH locus were reported to affect milk productions in cattle (Dybus et al. Reference Dybus, Grzesiak, Szatkowska and Błaszczyk2004), growth rate (An et al. Reference An, Hou, Wang, Li, Wang, Song, Zhou, Han, Ling and Cao2010) and milk yield (Malveiro et al. Reference Malveiro, Pereira, Marques, Santos, Belo, Renaville and Cravador2001; Dettori et al. Reference Dettori, Rocchigiani, Luridiana, Mura, Carcangiu, Pazzola and Vacca2013) in goats, milk yield in Serra da Estrela (Marques et al. Reference Marques, Santos, Carolino, Belo, Renaville and Cravador2006) and Sarda ewes (Vacca et al. Reference Vacca, Dettori, Balia, Luridiana, Mura, Carcangiu and Pazzola2013). Sequences of the complete GH1 (GenBank accession number: DQ450146-DQ450147), GH2-N (GenBank: DQ461644-DQ461681) and GH2-Z (GenBank: DQ461615-DQ461643) gene copies have been reported for the Serra da Estrela sheep, revealing several polymorphic sites (Marques et al. Reference Marques, Santos, Carolino, Belo, Renaville and Cravador2006; Marques, Reference Marques2007).
In this research we separately amplified, by PCR, the two copies of the sheep duplicate GH gene, which differ only in terms of their 5′ flanking regions. Then each gene copy was completely sequenced in all the sheep of the experiment. Analysis of the sequencing traces revealed polymorphisms occurring in each copy, and the relationship between polymorphisms and milk traits was evaluated.
The Sarda sheep breed is autochthonous to the island of Sardinia, where about 3 million head are reared and it is the most important Italian dairy sheep breed (Vacca et al. Reference Vacca, Carcangiu, Dettori, Pazzola, Mura, Luridiana and Tilloca2008). Owing to its high production and rusticity and adaptability to Mediterranean climate characteristics (Carcangiu et al. Reference Carcangiu, Luridiana, Vacca, Daga and Mura2011), it has spread to other Italian regions and to several North African Countries, where it is reared either as purebred or crossbred (Vacca et al. Reference Vacca, Dhaouadi, Rekik, Carcangiu, Pazzola and Dettori2010).
In the present investigation, the GH1/GH2-N and GH2-Z gene copies were entirely sequenced in a population of Sarda breed sheep (Ovis aries), in order to highlight sequence polymorphisms and evaluate association between genetic variants at the GH1/GH2-N and GH2-Z gene copies and milk traits.
Materials and methods
Animals and samples
One-hundred-and-six clinically healthy multiparous ewes of Sarda breed, aged 4 years, in their third lactation, were randomly selected from four commercial farms showing similar management and feeding systems, located on the island of Sardinia (Italy). Daily milk yield was recorded from each sheep of the experiment once a month, from February to June which corresponds to different stages of the milking period after lambs’ weaning according to ICAR guidelines (ICAR, 2012), measured as days in milk (DIM) from 15 ± 3 to 135 ± 3 d (15, 45, 75, 105 and 135 DIM). Individual milk samples were taken for milk analysis and individual blood samples were taken for DNA extraction.
Milk was collected by hand milking. Approximately 5 ml of blood was collected by puncture via the jugular vein with only physical restraint. Blood samples were taken solely for research purposes and experienced veterinarians performed blood collection and all efforts were made to minimise suffering. No specific permission to an animal ethics committee was required for these locations and activities, because according to Article 2 of the EC Directive 86/609/EEC, none of the procedures used here met the criteria to be defined as experiments, then specific approval was not necessary.
Milk analysis
Milk traits recorded included milk yield; protein, fat and lactose percentage, determined by infrared spectrophotometer (Milko-Scan 133B; Foss Electric, DK-3400 Hillerød, Denmark) (IDF Standard No 141C, 2000); casein content (IDF Standard No 29, 1964); milk energy, estimated on the bases of calories provided by total solid components (National Research Council, 2007).
DNA analysis
DNA was obtained from leucocytes using a commercial kit (NucleoSpin Blood, Macherey-Nagel). DNA concentration and purity were measured by spectrophotometer (Eppendorf Biophotometer, Hamburg, Germany). The GH1/GH2-N gene was amplified by polymerase chain reaction (PCR) with the primer pair GHTF/GHTR (Marques et al. Reference Marques, Santos, Carolino, Belo, Renaville and Cravador2006) to produce a DNA fragment of about 2·05 kb, and the GH2-Z gene was amplified with the primer pair GH2ZF/GH2ZR to give a fragment of about 2·22 kb (Vacca et al. Reference Vacca, Dettori, Balia, Luridiana, Mura, Carcangiu and Pazzola2013) (Table 1). The two gene copies were analysed by Sanger sequencing in all the animals of the experiment, in the forward direction. Sequencing oligonucleotides were designed with the Primer3Plus software (http://www.bioinformatics.nl/cgi-bin/primer3plus/primer3plus.cgi/) based on GenBank: DQ450546 and DQ461643 (Table 1). Prior to sequencing, amplicons were purified with Agencourt® AMPure® kit (Beckman Coulter, USA), then sequencing was obtained with an Applied Biosystems 3730 DNA Analyser (Applied Biosystems, Foster City CA, USA). To test the quality of sequencing reactions, about 15% of the samples were also sequenced in the reverse direction, as well as those subjects who showed a heterozygous insertion/deletion (indel), and thus alteration of the sequencing reaction downstream of it.
Oligonucleotides utilised for PCR amplification and sequencing, at the sheep GH locus
Analysis of the sequencing data
The Phred/Phrap/Crossmatch softwares were used to compute data quality information from chromatogram files and to construct contig sequences (Ewing & Green, Reference Ewing and Green1998; Ewing et al. Reference Ewing, Hillier, Wendl and Green1998). PolyPhred was used to compare sequence traces and search for homozygotes and heterozygotes in these contigs (Nickerson et al. Reference Nickerson, Tobe and Taylor1997; Bhangale et al. Reference Bhangale, Stephens and Nickerson2006). PolyPhred polymorphism tags were examined using Consed for each putative SNP (Gordon et al. Reference Gordon, Abajian and Green1998; Reference Gordon, Desmarais and Green2001; Gordon, Reference Gordon, Baxevanis and Davison2004). Sequence polymorphisms were described according to the Human Genome Variation Society (http://www.hgvs.org/mutnomen/).
The Haploview software package (Barrett et al. Reference Barrett, Fry, Maller and Daly2005) was used to estimate minor allele frequency (MAF), observed and expected heterozygosity, and to test genotype distributions for departure from Hardy-Weinberg equilibrium at each polymorphic locus.
The Alibaba (v.2·1) software and the TRANSFAC (v.7·0) database (http://www.gene-regulation.com/pub/programs/alibaba2) were used to find motifs in promoter regions. Nucleotide sequences were translated into amino acids using ExPASy Bioinformatics Resource Portal (http://www.expasy.org/) and compared with the translated amino acid sequence of GenBank: DQ450146 (GH1) and GenBank: DQ461643 (GH2-Z).
Association analysis
Association analysis between the SNPs at GH1/GH2-N and GH2-Z and milk traits (milk yield, protein, fat, casein and lactose percentages and energy) was performed with the SAS Software (SAS Inst. Inc., Cary NC, USA) using the mixed effect model procedure and the following model:
$${y_{ijkl}} = {\rm \mu} + {H_i} + {\rm DI}{{\rm M}_j} + {G_k} + H \times {\rm DI}{{\rm M}_{ij}} + {a_{\rm l}}({G_k}) + {e_{ijkl}}$$where y ijkl is the analysed variable, μ is the general mean, H i is the fixed effect of the herd (i = 4 levels), DIMj is the fixed effect of the stage of milking measured as days after lambs’ weaning (j = 5 levels 15, 45, 75, 105 and 135 DIM), G k is the fixed effect of genotype at each SNP (k = 2 or 3 levels), H × DIMij is the interaction effect between the herd and the stage of lactation, a l(G i) is the random effect of lth animal nested within the kth genotype and e ijkl is the error effect. Because of the low frequency of some genotypes and the resulting rank deficiency in contingency tables, each SNP was analysed at a time. In order to avoid a bias of results, SNPs with MAF <0·10 were excluded from the model. Bonferroni method was used to correct multiple comparison of the means and significance was set at P < 0·05.
Results
Analysis of the sequencing data
PCR amplification of the GH2-Z gene was positive in all the ewes, which means that no subject carried the Gh1 allele in homozygosis; however, as it cannot be excluded the occurrence of Gh1/Gh2 heterozygotes, the GH2-N gene was always considered together with GH1 (GH1/GH2-N). Alignment of the sequencing traces from 106 ewes allowed the identification of 72 nt changes, 39 in the GH1/GH2-N gene and 33 in the GH2-Z gene (Tables 2 and 3). All the nt changes, as well as population genotypes, were submitted to the dbSNP database (https://www.ncbi.nlm.nih.gov/SNP/) and were given accession numbers: ss748770506-ss748770580.
SNPs at the GH1/GH2-N gene in Sarda sheep
Polymorphic sites shared by the GH1/GH2-N and GH2-Z gene copies are shown in bold type. SNP ID, submitted SNP accession number from the dbSNP database; DQ450146, nucleotide positions relative to GenBank DQ450146 (GH1); X12546, nucleotide positions related to GenBank X12546; MAF, minor allele frequency; Ho, observed heterozygosity; He, expected heterozygosity; HWpval, Hardy Weinberg P value
† Variations already described in Vacca et al. (Reference Vacca, Dettori, Balia, Luridiana, Mura, Carcangiu and Pazzola2013)
‡ Variations already described in Marques et al. (Reference Marques, Santos, Carolino, Belo, Renaville and Cravador2006)
§ Variations already described in Marques (Reference Marques2007)
¶ Single nucleotide deletion
†† Rare allele in brackets
SNPs at the GH2-Z gene in Sarda sheep
Polymorphic sites shared by the GH1/GH2-N and GH2-Z gene copies are reported in bold type. SNP ID, submitted SNP accession number from the dbSNP database; DQ461643, nucleotide positions relative to GenBank DQ461643 (GH2-Z); X12546, nucleotide positions in GenBank X12546; MAF, minor allele frequency; Ho, observed heterozygosity; He, expected heterozygosity; HWpval, Hardy Weinberg P value
† Variations already described in Vacca et al. (Reference Vacca, Dettori, Balia, Luridiana, Mura, Carcangiu and Pazzola2013)
‡ Variations already described in Marques et al. (Reference Marques, Santos, Carolino, Belo, Renaville and Cravador2006)
§ Variations already described in Marques (Reference Marques2007)
¶ Single nucleotide deletion
†† Rare allele in brackets
At the GH1/GH2-N gene, 27 nt changes were transitions, 7 transversions and 5 single nucleotide indels, while at the GH2-Z gene, 24 nt changes were transitions, 7 transversions and 2 were indels. A total of 11 nt changes were detected within coding regions: 3 occurred at the GH1/GH2-N exons 2 and 3 (therefore exons 1, 4 and 5 were monomorphic) and 8 coding SNPs occurred at the GH2-Z gene (only exon-1 was monomorphic). All the coding SNPs derived from base transitions (Table 4). Alignment of all the GH1/GH2-N sequencing traces to GenBank DQ450146, representing the complete GH1 gene sequence, revealed two further nt changes: ss748770542 and ss748770543, which were monomorphic in the population study. Alignment of the GH2-Z sequencing traces to GenBank DQ461643, which is the complete sequence of the sheep GH2-Z gene, revealed one further nt change ss748770569, monomorphic in the population study. Genotype distributions were in Hardy-Weinberg (HW) equilibrium at all loci at the GH1/GH2-N gene, while 13 SNPs did not respect HW equilibrium at the GH2-Z gene.
Coding SNPs and putative amino acid changes at the Sarda sheep GH gene
Association analysis
Association was evaluated between all the SNPs with MAF higher than 0·10 and milk traits (milk yield, fat, protein, casein and lactose percentages and energy) and values were reported in Tables 5 and 6. No interaction effects on milk traits were observed for DIM and genotype. Statistical analysis revealed a significant effect of herd, as well as a significant effect of the stage of lactation, on all milk traits under study.
Analysis of variance (F values) showing the effects of herd (H), days in milking (DIM) and genotype (G) on milk traits
*, **, and *** indicate significant F-values at P < 0·05, 0·01 and 0·001, respectively; ns not significant
LS means of milk traits, according to genotypes at sheep GH gene
n, number of animals.
* Different capital letters in the same column indicate values significantly different (after Bonferroni correction) for each SNP at P < 0·01; lower case letters indicate significance at P < 0·05
Significant differences in milk traits were observed among genotypes at polymorphic loci only at the GH2-Z gene (Tables 5 and 6). Sheep with homozygote genotype ss748770547 CC had higher fat percentage and energy than TT (P < 0·01); ss748770547 TT showed the lowest casein (P < 0·05) and lactose (P < 0·01) percentage. Sheep with genotype ss748770553 GG had higher milk yield (P < 0·05) and lactose percentage (P < 0·01) than AG and AA, while ss748770553 AA sheep showed lower milk protein and casein percentages (P < 0·01) and lower energy values (P < 0·05) than AG and GG. Homozygotes ss748770572 GG sheep had higher milk yield and lactose percentage than GT (P < 0·01) and was associated with lower fat percentage than GT and TT (P < 0·01).
Discussion
Growth hormone, encoded by the GH gene, modulates growth and lactation in domestic animals, and then genetic variants at the GH locus might affect its level, stage or tissue of expression, as well as its functions, having an impact on production traits. No homozygote Gh1 subjects were found in the study population, while 1·14% Gh1/Gh1 subjects have been reported for the Serra da Estrela sheep (Marques et al. Reference Marques, Santos, Carolino, Belo, Renaville and Cravador2006). Also the Texel sheep, utilised for the sheep genome reference sequence, carries the Gh1 allele, mapped to chromosome 11, GenBank: NC_019468·1 (47540064.47541799, complement) (International Sheep Genomics Consortium et al. Reference Archibald, Cockett, Dalrymple, Faraut, Kijas, Maddox, McEwan, Hutton Oddy, Raadsma, Wade, Wang, Wang and Xun2010). The Gh1 allele has been found to be associated with lower milk yield in the Serra da Estrela sheep (Marques et al. Reference Marques, Santos, Carolino, Belo, Renaville and Cravador2006), then this allele might have undergone negative selection in the Sarda sheep, which has been selected for many decades to improve milk production (Mura et al. Reference Mura, Daga, Paludo, Luridiana, Pazzola, Bodano, Dettori, Vacca and Carcangiu2012).
A total of 72 polymorphic sites were evidenced at the GH1/GH2-N and GH2-Z genes, plus 3 monomorphic nt changes compared to GenBank: DQ450146 and DQ461643 (dbSNP ss748770506 to ss748770580). The two sheep GH gene copies shared seven common SNPs, among these, the coding SNP c.255G>A (ss748770530 and ss748770563) and the intronic SNP c.453 + 35A>G (ss748770534 and ss748770570) have been reported also at the caprine GH locus (Dettori et al. Reference Dettori, Rocchigiani, Luridiana, Mura, Carcangiu, Pazzola and Vacca2013), indicating that they probably occurred in a common ancestor. Also SNP ss748770564, occurring only at sheep GH2-Z, has been revealed in the caprine GH gene, indicating that some sequences of the GH gene are more susceptible than others to the occurrence of mutations. According to MacDonald et al. (Reference McDonald, Wang, Huang and Leu2011), indels occur in DNA regions with polynucleotide repeats or other repeat sequences, which are associated with the accumulation of nucleotide substitutions, over evolutionary time scales. For example, the indel c.*30delT (sheep ss748770538 and ss748770577) also occurs in goats, but with an inserted/deleted C (c.*30delC) rather than T (Vacca et al. Reference Vacca, Dettori, Balia, Luridiana, Mura, Carcangiu and Pazzola2013) and we detected a total of five indels at GH1/GH2-N and two at GH2-Z, similar to the ones revealed in the caprine GH gene (Dettori et al. Reference Dettori, Rocchigiani, Luridiana, Mura, Carcangiu, Pazzola and Vacca2013).
Some authors report that amino acid positions 9 and 63 (corresponding to 35 and 89 of the entire protein sequence, respectively) of the sheep GH protein sequence, allow differentiating the two gene copies GH2-N and GH2-Z (Lacroix et al. Reference Lacroix, Devinoy, Servely, Puissant and Kann1996; Wallis et al. Reference Wallis, Lioupis and Wallis1998). We found that GH1/GH2-N always showed a glycine at position 35, while GH2-Z had an arginine, then this amino acid position was discriminating. In contrast, amino acid position 89 was not discriminating, as GH2-Z was polymorphic (p.S89G).
The two copies of the sheep GH gene differ essentially in their 5′ flanking sequences (Marques et al. Reference Marques, Santos, Carolino, Belo, Renaville and Cravador2006), then we analysed these sequences to check whether the nt changes therein were involved in the alteration of transcription factor binding sites (TFBS). Several SNPs were part of putative TFBS, among them, SNP ss748770506 was part of a potential binding site for the pituitary-specific transcription factor POU1FI (formerly Pit-1), which remained unchanged when G or A occurred at nt position c.−183. The POU1FI transcription factor regulates the expression of pituitary genes, such as growth hormone and prolactin genes (Mura et al. Reference Mura, Daga, Paludo, Luridiana, Pazzola, Bodano, Dettori, Vacca and Carcangiu2012), specifically binding the consensus sequence 5′-TAAAT-3′. Two sequences potentially recognised by this TF were found at the GH1/GH2-N gene copy (−183/−174 and −142/−133), but no such sequence was found at the promoter of GH2-Z, which probably, during the duplication process, has retained only 82 nt of the 5′ flanking region of the original GH1 gene, losing the two potential POU1FI binding sites. This might explain the different tissues of expression of the two genes, pituitary for GH1/GH2-N and placenta for GH2-Z (Lacroix et al. Reference Lacroix, Devinoy, Servely, Puissant and Kann1996). A thorough analysis of the promoter regions at the sheep GH locus has been studied by Marques (Reference Marques2007).
Association analysis
Association analysis showed that polymorphism of three SNPs at the GH2-Z gene affected milk traits. SNP ss748770547, occurring exclusively in the 5′ flanking region of GH2-Z (c.−259 T>C), significantly affected milk fat percentage, as well as lactose percentage and energy values. This SNP was part of a potential transcription factor binding site (TFBS) for C/EBP alpha, when the C occurred, but when T substituted for C this binding site was lost, as already evidenced in the same sheep breed (Vacca et al. Reference Vacca, Dettori, Balia, Luridiana, Mura, Carcangiu and Pazzola2013). Transcription factor C/EBP alpha (CCAAT/Enhancer Binding Protein) is known to be involved in the regulation of adipogenesis and adipoblast differentiation (Payne et al. Reference Payne, Au, Lowe, Rahman, Friedman, O'rahilly and Rochford2009; Siersbaek et al. Reference Siersbaek, Nielsen and Mandrup2012), then ss748770547 may be a causal mutation affecting milk fat content, which has been highlighted here for the first time.
It is known that pituitary GH, in ruminants, is involved in the mobilisation of body lipid stores, in the early weeks of lactation, by an increase of lipoprotein lipase activity in the adipose tissues, with an increase of free fatty acids in plasma, which become available to the mammary gland (Palmquist, Reference Palmquist, Fox and McSweeney2006; Vijayakumar et al. Reference Vijayakumar, Novosyadlyy, Wu, Yakar and LeRoith2010). In addition, pituitary GH can affect the synthesis of fatty acids de novo via activation of Acetyl-CoA carboxylase, mediated by IGF-1 (Palmquist, Reference Palmquist, Fox and McSweeney2006). The mechanism through which the sheep placental GH, expressed from day 27 to day 55 of lactation, may be involved in a process of adipogenesis, remains to be elucidated, and it could be different from that of pituitary GH.
Polymorphism at the coding SNP ss748770553 affected milk yield, as well as protein, casein and lactose contents. This SNP is located in exon-2 of GH2-Z and corresponds to the nonsynonymous change p.Ala27Thr, which is the N-terminal residue of the mature protein. Although specific studies on the regulation of the half-life of GH are not reported in the literature, the abovementioned N-terminal residue might be involved in the regulation of GH degradation, through the N-End Rule Pathway (Tasaki et al. Reference Tasaki, Sriram, Park and Kwon2012).
SNP ss748770572 (intron 4 of GH2-Z) significantly affected milk yield and lactose content (P < 0·01). It should be noted that SNP ss748770572, describing the nt change c.453 + 84 T>G, was common to both GH1/GH2-N (SNP ss748770535) and GH2-Z, but ss748770535 was excluded from the statistical analysis because it showed MAF <0·10. The effect of this intronic SNP might be due to a modification of the splicing mechanisms within the GH locus. In fact, alternative splicing has been observed in the human GH gene cluster, including five genes which, despite the high degree of sequence identity, show differences in tissue distribution, expression levels, but also in patterns of splice site selection (Männik et al. Reference Männik, Vaas, Rull, Teesalu and Laan2012).
Conclusions
GH plays an important role in the metabolism of mammals and can affect livestock production. In this study we revealed that the duplicate sheep GH gene shows high polymorphism. The promoter region of GH1/GH2-N, expressed in the pituitary, contains two potential binding sites for the pituitary-specific transcription factor POU1FI, which are lost at the GH2-Z gene, expressed only in the placenta. Statistical analysis revealed that polymorphism of SNPs at the GH2-Z gene affected milk traits. SNP ss748770547 significantly affected fat percentage. This information might be useful in selection schemes for dairy sheep breeding.
We gratefully acknowledge the Regione Autonoma della Sardegna (L.R.7/2007 CUP: J71J10000290002 and L.R.7/2007: MIGLIOVIGENSAR) for financial support.
