Associations of myostatin gene polymorphisms with performance and mortality traits in broiler chickens

Myostatin is a negative regulator of skeletal muscle growth. We evaluated effects of myostatin polymorphisms in three elite commercial broiler chicken lines on mortality, growth, feed conversion efficiency, ultrasound breast depth, breast percentage, eviscerated carcass weight, leg defects, blood oxygen level, and hen antibody titer to infectious bursal disease virus vaccine. Progeny mean data adjusted for fixed and mate effects and DNA from 100 sires per line were used. Single nucleotide polymorphisms (SNPs) of the myostatin gene segregating in these lines were identified by designing specific primers, amplifying individual DNA in each line by polymerase chain reaction, cloning, sequencing and aligning the corresponding products. Individual sires were genotyped for five identified SNPs which contributed to eight haplotypes. Frequencies of SNP alleles and haplotypes differed between lines. Using the allele substitution effect model, the myostatin SNPs were found to have significant (P < 0.031) associations with growth, mortality, blood oxygen and hen antibody titer to infectious bursal disease virus vaccine, although the associations were not often consistent across lines. These results suggest that the myostatin gene has pleiotropic effects on broiler performance.


INTRODUCTION
Myostatin is a member of the transforming growth factor-beta (TGF-beta) superfamily and highly conserved in gene structure among vertebrate species [15,16]. The myostatin gene is highly polymorphic. Nineteen single nucleotide polymorphisms (SNPs) and 20 haplotypes have been detected in 28 European bovine breeds, of which some were breed-specific [6]. Baron et al. [1] identified seven SNPs and one deletion in exon 2 of the myostatin gene in broiler and/or layer chicken lines. Gu et al. [12] identified three SNPs in the 5 -regulatory region and two SNPs in the 3 -regulatory region in the chicken, and these differed in allele frequencies between breeds.
Myostatin is mainly expressed in muscular tissues. Several studies have reported that myostatin negatively regulates skeletal muscle growth. Muscle weights of mice lacking the myostatin gene are about two to three times greater than those of wild type mice [18]. A loss-of-function mutation in the myostatin gene causes double muscling phenotypes in Belgian Blue and Piedmontese cattle [17]. Myostatin negatively regulates activation of satellite cells that promote postnatal muscle growth and repair in mice [16]. Overexpression of myostatin during infection with chronic infectious respiratory diseases has been shown to reduce skeletal muscle growth in pigs [7]. In humans, a loss-of-function mutation in the myostatin gene in children increases muscle bulk and strength [21]. The function of myostatin in negative regulation of skeletal muscle growth is modulated by additional genetic factors, as illustrated by the finding that myostatin protein levels are not significantly different between broilers and layers [11]. Other genes have been shown to affect muscularity in the presence of Mstn Cmpt-dl1Abc (a 12-bp deletion in the myostatin gene) in mice [23]. Potts et al. [19] reported 19 genes to be differentially expressed between normal and double-muscled bovine embryos.
Myostatin gene polymorphisms have effects on multiple traits. Mutations in myostatin regulatory regions have been shown to be associated with abdominal fat weight, abdominal fat percentage, birth weight, breast muscle percentage and breast muscle weight in chickens [12]. Furthermore, Casas et al. [5] found that calves with two copies of the inactive allele of myostatin were more likely to die before weaning, were heavier at birth and leaner, and had more muscle mass than animals with zero or one copy; heterozygous calves had the heaviest bodyweight at weaning and the highest live weight; calves with zero copy were the highest in fat content.
The objectives of the current research were to identify genetic polymorphisms in the myostatin gene in elite commercial broiler chicken lines and to determine their associations with growth, mortality, yield and support traits in high-and low-health environments. Results of this study will increase our understanding of the many functions of the myostatin gene.

Populations and traits
About 100 sires from each of three elite commercial broiler chicken lines (coded X, Y and Z) were used. Mortality and performance data were recorded on 17 to 437 progeny from each sire (Tab. I). A sire's progeny were raised in one high hygiene environment (HH) or in one of two low hygiene environments (LH: LH1 and LH2). Traits recorded (Tab. I) in all three environments were body weight at 7 and 40 days of age (BW7 and BW40) and mortality from hatching to 14 days (EMORT) and from 14 to 40 days of age (LMORT). Nine additional traits were recorded in the HH environment: ultrasound breast depth (US), percentage of breast (BR), feed conversion efficiency (FCR), eviscerated carcass weight (EV), three leg-defect related traits (measured on a subjective scale by trained technicians) [twisted legs or evident tibial dischondroplasya (Leg), X-ray-inspection-based sub-clinical or incipient development of tibial dischondroplasya (Lixi), and curly or crooked toes and/or bowed legs (Tobo)], and blood oxygen content measured by pulse oximeter (Oxi) and female's antibody titer to infectious bursal disease at 27 wks (IBD).

DNA samples, identification of SNPs and genotyping
Blood samples from each sire were individually collected in tubes containing EDTA and immediately frozen at −20 • C. Individual genomic DNA was extracted from a 15 µL aliquot of blood with DNAzol reagent (Invitrogen, Carlsbad, California, USA) following the manufacturer's protocol. DNA samples were quantified by spectrophotometer and visualized in agarose gel. One hundred ng of genomic DNA was used as a PCR template.
Five pairs of gene-specific primers were designed based upon published myostatin gene sequence (GenBank accession number: AF346599). Primer sequences and their annealing temperature are in Table II and were used to amplify the coding region, part of the promoter and intron regions of myostatin. Three primer pairs, MSTpr, MSTex1 and MSTex3, were used to identify and genotype SNPs. Primers MSTex2-1 and MSTex2-R were used to identify the SNP that was genotyped using primers MSTex2-2 and MSTex2-R. Independent amplifications were conducted for individual DNA samples from each line using a high fidelity Taq DNA polymerase (Invitrogen, Carlsbad, California, USA). PCR products were cloned using Topo Cloning Kit (Invitrogen, Carlsbad, California, USA) and 192 plasmids clones from each line and myostatin region were sequenced using the BigDye terminator kit on an Table I. Phred/Phrap and Consed programs were used for sequence quality evaluation, contig assembling and SNP inspection, as previously described by Gordon et al. [8]. Clustal-W software [22] and Vector NTI (Invitrogen, Carlsbad, California, USA) were used for sequence alignments, translation and polymorphism identification. Individual sires were genotyped for five of the identified SNPs using PCR-RFLP methods. These five SNPs were chosen based on the feasibility and cost of genotyping using PCR-RFLP, partially in view of their potential future routine application in breeding programs. The PCR reaction was performed in 55 µL of reaction solution with 0.025 units/µL Taq (Invitrogen Ltd., Paisley, UK), 1.5 mM Mg, 200 µM dNTPs, 0.5 µM primers, and 5 µL DNA sample (25 ng/µL). Cycling conditions for PCR were 94 • C 2 min for 1 cycle, 94 • C 30 s, annealing temperature (Tab. II) 30 s and extension 72 • C 30 s for 40 cycles, and extra extension 72 • C for 2 min. Digestion of PCR products was performed in a 12 µL of reaction solution following manufacturer's directions (New England Biolabs, Ipswich, MA, USA) and other conditions (Tab. II). DNA fragments were electrophoretically separated on 3% Metaphor gel (Cambrex Bio Science Wokingham, Ltd., Berkshire, UK).

Statistical analysis
Associations of polymorphisms with traits were identified based on analysis of the mean adjusted performance of progeny of each sire in a given environment. These were derived from progeny performance records that were each adjusted for the systematic environmental effects of sex, age, hatch and mating group, and for mate effects (half the EBV of the dam) by subtracting solutions obtained from routine animal model genetic evaluation procedures that are used within these lines. Association analyses were conducted using weighted least squares of SAS  PROC GLM (SAS  9.1), with weights equal to the number of progeny included in the adjusted mean to account for differences in variance of residuals. Only additive associations were evaluated because progeny means primarily reflect the additive effects of genes, being related to half the sire's breeding value. Association analyses were performed separately for single SNPs and for haplotypes of SNPs. In addition, an outlier analysis was performed. Twenty-one possible outlier data points were identified in the 14 traits by visual inspection of distributions of progeny means against numbers of progeny. Comparison of results from analyses with all data points with those after removal of the outliers showed that the identified outliers had minor  effects on the SNP-trait association results. The results reported in this paper are, therefore, based on all data points in the original dataset. In all analyses, residuals were assumed uncorrelated. Genetic relationships among sires were limited and were ignored in the analyses; the 102, 100, 114 sires that were analyzed for lines X, Y, and Z originated from a total of 51, 50, and 66 sires and 71, 75, and 93 dams. Thus, half-and full-sib relationships among the evaluated sires were small and ignoring them is not expected to bias results.
Analysis with single SNPs: Three allele substitution effect models were employed to determine associations between an SNP and a trait. Model 1 was for within-line data analysis and Models 2 and 3 were for analyses across lines. where, y i and y i j are the adjusted progeny mean of the ith sire of line j; µ is a general mean; L j is the effect of the jth line (j = 1, 2, 3); f i and f ji represent the number copies of a given allele of the SNP carried by the ith sire; b j is the allele substitution effect for line j; ε i and ε ji are the residual for sire i with variance = 1 N i σ 2 e (N i is the number of progeny for sire i). Significance level for addressing an association was set to be 0.031 corresponding to a 50% proportion of false positive associations (PFP).

Analysis with haplotypes:
It is likely that the SNPs within the myostatin gene are in linkage disequilibrium because of their proximity on the chromosome, and the haplotypes formed by these SNPs may be associated with traits of interest. Thus, presence of linkage disequilibrium among myostatin SNPs was determined using a χ 2 test, amount of linkage disequilibrium was measured using r 2 [13], and haplotype frequencies were estimated by maximum likelihood using the software Arlequin version 2.000 [20]. For sires whose haplotypes could not be inferred with certainty, haplotype frequencies were used to assign haplotype probabilities. Models to analyze associations of haplotypes with traits were the same as described above but replacing SNP allele effects by n−1 h=1 b jh f i jh , where f ih is the probability of haplotype h in sire i (or the sum of the probabilities that the first and second haplotype of the sire are h) and b h is the substitution effect for haplotype h. In these models, the effect of the most frequent haplotype (GAGCT for each line) was set equal to zero, such that the models were of full rank. As a result, b h represents the effect of substituting a copy of haplotype GAGCT by a copy of haplotype h. Significance tests for haplotype effects on traits were performed in two steps: first, the overall haplotype effect for a trait was tested; if this was significant, then the contrast of each haplotype with the most frequent haplotype was tested.

SNPs and haplotypes
Thirteen SNPs were identified in the amplified regions of exons 1, 2 and 3 and introns 1 and 2 in myostatin gene (Tab. III). Of them, five exonic SNPs (MST2100, MST2109, MST2244, MST2373 and MST4842) were genotyped on individual sires using the PCR-RFLP method to determine allele frequencies and associations with phenotypic traits. Allele frequencies of the five SNPs varied by line (Tab. IV). The MST4842 SNP was fixed in lines Y and Z. All SNPs were in Hardy-Weinberg equilibrium (P > 0.1) but in linkage disequilibrium with each other (P < 0.01) (Tab. IV). The SNPs at positions 2100, 2109, and 2373 bp in the myostatin gene sequence (GenBank accession number: AF346599) were in complete disequilibrium in lines Y and Z, but an intervening SNP at 2244 bp was in incomplete disequilibrium (r 2 < 0.84) with the other SNPs. A total of eight haplotypes segregated in line X but only three of these segregated in lines Y and Z (Tab. IV), which was much fewer than the 32 haplotypes that were possible. Haplotypes GAGCT and GGCCG were prevalent in line X and haplotypes GAGCT and AGCTT were prevalent in lines Y and Z (Tab. IV). Haplotype GAGCT was the most frequent haplotype in all three lines.

Trait associations of single SNPs
Results for associations of individual myostatin SNPs with traits are in Figure 1. Results in Figure 1 were obtained from analysis with Model 1. Because MST2373 was in complete disequilibrium with MST2109 in lines Y and Z, results for the SNP are not shown. Compared with results from Model 1, analyses with Models 2 and 3 seldom provided extra information and are, therefore, not shown. For analyses with Model 2, significant interactions between SNPs and lines were detected when allele effects were significant in at least one line in Model 1. In these cases, the results for individual lines from Model 1 are reported. Analysis with Model 3 detected only one new association, for MST2100 with Tobo. MST2100: Allele 1 (nucleotide A) of MST2100 had consistent negative effects on growth in the HH environment across ages and lines, although its effect was significant only on BW7 in line Z (Fig. 1A). Allele 1 also had a negative effect on IBD in lines X and Z and was significant (P < 0.031) in line Z. Allele 1 had a non-significant positive effect on Oxi in line Z but a significant (P < 0.001) negative effect on Oxi in line X.

MST2109:
In lines Y and Z, the effects of allele 1 on the traits were similar in magnitude to those for MST2100 because of the complete LD between these two SNPs, but the direction of effects was opposite (Fig. 1B). In line X, allele 1 had consistent positive effects on growth in the LH environments, with two of four traits having significant effects (P < 0.031). Allele 1 also had a negative effect on Leg (P < 0.04), a positive effect on Oxi (P < 0.043), and positive effects on LMORT in the LH environments (Fig. 1B).

MST2244
: Allele 1 had consistent negative effects on growth in the LH environments in line X and Z and was significant for three of the four growth traits in line X (P < 0.031) (Fig. 1C). Allele 1 had a consistent negative effect on LMORT, with one allele effect being significant (P < 0.01) in the LH environments in line X (Fig. 1C).
MST2373: MST2373 had the same associations with the traits as MST2109 in lines Y and Z because of complete linkage disequilibrium (Tab. IV). In line X, allele 1 had a positive effect on mortality in the HH environment, a positive effect on IBD (not significant), and a significant (P < 0.001) positive effect on Oxi (Fig. 1D).

MST4842:
Associations of MST4842 with traits were evaluated only in line X because of its fixation in lines Y and Z. In line X, allele 1 had consistent negative, though partially significant (P < 0.031), effects on growth across ages in the LH environments (Fig. 1D). Allele 1 also had negative effects on LMORT in LH environments, though only significant (P < 0.01) in LH2.

Associations of haplotypes with traits
Few significant (P < 0.1) haplotype-trait associations were detected (Tab. V) and they generally matched results from the individual SNP analyses. In line X, haplotypes of the five SNPs had consistent associations (P < 0.029) with LMORT across environments and ages. Haplotype GAGCT had significantly lower LMORT than haplotype AGGTG in both LH and HH environments. Haplotype GAGCT also had significantly lower EV than AGCTT and GAGCG, but significantly higher EV than AGGTT. In line Y, there were no significant associations of haplotypes with the traits studied. In line Z, myostatin haplotypes had significant associations with BW7-LH1, BW7-HH, IBD-HH, and LMORT-LH2 (Tab. V). Of the three haplotypes in line Z, AGGTT had the highest BW7 in the LH1 environment but the lowest BW7 in the HH environment; GAGCT had the highest IBD-HH; and AGCTT had the lowest LMORT in the LH2 environments.

DISCUSSION
The high degree of polymorphism of the DNA sequence of the myostatin gene, which has been demonstrated in previous studies [1,6,12], was further displayed here in three elite commercial broiler chicken lines. Of the identified 13 SNPs, 11 were in exons and 2 were in introns. Four SNPs, MST2283, MST2346, MST2416 and MST7434, resulted in disappearance/appearance of restriction sites that were expensive for PCR-RFLP genotyping. Four SNPs, MST4405, MST4954, MST7435 and MST7436, did not result in disappearance/appearance of restriction sites. The remaining five of the 13 SNPs were genotyped using PCR-RFLP. The four genotyped SNPs, MST2100, MST2109, MST2244 and MST2373 in exon 1 were synonymous substitutions (GCG → GCA, CCG → CCA, CGC → CGG, GAC → GAT and ACG → ACT, respectively). Therefore, these four SNPs are genetic markers and likely not causative mutations for the traits with which they were found to be associated. The SNP MST2373 was the same as that previously reported (as nt. 1204 C/T) by Baron et al. [1] in two chicken lines that were separately selected for bodyweight and egg production [1]. The SNP MST4842 in exon 2 was a non-synonymous substitution (CTG → CGG), which resulted in a change of an amino acid (Leu → Arg). Baron et al. [1] did not identify MST4842 but found six other non-synonymous SNPs in exon 2. The five SNPs found in the current study formed eight haplotypes in the broiler chickens; line X had five line-specific haplotypes, of which haplotype GGCCG was most frequent.
Each of the five identified myostatin SNPs was significantly (P < 0.031) or consistently associated with bodyweight at day 7 and/or day 40 in at least one broiler line. Associations of the non-synonymous SNP MST4842 with bodyweight in line X may be due to the amino-acid change in myostatin, although functional studies will be needed to ascertain this. Associations of myostatin polymorphisms with growth have been frequently reported in humans and cattle [2-5, 9, 10, 14, 15, 17, 21]. These imply that, because the main function of myostatin is regulation of skeletal muscle growth, changes in its protein sequence can modify growth function. The detected associations of the synonymous SNPs with growth may be attributed to their linkage disequilibrium with other polymorphisms that exist in regulatory regions or exons of myostatin or to QTL beyond the myostatin gene. Gu et al. [12] reported a polymorphic site in the 5 regulatory region of myostatin, which was associated with hatch weight in F 2 chickens from a broiler by Silky cross [12]. This or other regulatory-region polymorphisms could be in linkage disequilibrium with some of the synonymous SNPs detected in this study, which would explain their significant associations. Identification and further analysis of such polymorphisms were, however, beyond the scope of the present study.
The current study on myostatin SNP associations in chickens measured more extensive and diverse traits than have been reported in most previous studies. In addition to associations with growth, myostatin polymorphisms were also found to be associated with mortality, yield and support traits. This pleiotropic effect of the myostatin gene has also been reported by others in chicken and other species. Gu et al. [12] found homozygous genotypes AA and BB at a locus in the 5 regulatory region to be associated with higher abdominal fat weight and abdominal fat percentage than AB in the F 2 chickens from a cross of broiler and Silky chickens. Casas et al. [5] showed that calves with two copies of the inactive allele of myostatin were more likely to die before weaning, that calves with one copy of the inactive allele had leaner and more muscled carcasses than animals without inactive alleles, and that calves with zero copy of the inactive allele were highest in fat content.
Polymorphisms for myostatin have been used in marker-assisted selection in beef cattle but not in chickens. The current results suggest that one nonsynonymous SNP (MST484) and four synonymous SNPs may have pleiotropic associations with a set of traits and show genotype by environment interactions. This represents the first full evaluation of pleiotropic effects of myostatin mutations on different traits, but these results require confirmation in larger data sets and other lines. In addition, the causative reasons for the associations are not clear, even for the non-synonymous SNP MST484.
Pleiotropic effects and presence of genotype by environment interactions of myostatin can complicate genetic improvement but knowledge about them will allow breeding programs to be structured to minimize their impact or maximize their benefits. Presence of pleitropic effects will require identification of alleles or haplotypes that have a net beneficial effect across traits. Presence of genotype by environment interactions of myostatin implies that the network that the gene is involved in is sensitive to environmental changes. Knowledge of such genotype by environment interactions is important for genetic selection programs, to ensure that the allele(s) selected for has(have) beneficial effects in the target environment.