Sequence and polymorphism analysis of the camel ( Camelus dromedarius ) myostatin gene

Myostatin (MSTN), a negative regulator of skeletal muscle development in mammals, represents a key target for genetic investigations in meat-producing animals, with mutations responsible for increased skeletal-muscle mass currently described in several livestock species. Dromedary camels play a major economic role as suppliers of meat for human consumption across several countries. Notwithstanding, a comprehensive characterization of the sequence variability at the Camelus dromedarius MSTN locus was still lacking. Here we present the first extensive sequence and polymorphism analysis of the MSTN gene in the C. dromedarius species. Out of more than 3.6 kb of nucleotide sequence screened on 22 animals from 3 different Northern African regions, only 3 variant sites in the first intron were detected. The low observed diversity may reflect the evolutionary history of the species, likely developed as domesticates from a low variable wild ancestor population. Sequence identity among C. dromedarius and other Cetartiodactyla highlighted a tree topology consistent with previous reports of a closer relationship between Tylopoda and Suiformes. A close similarity between C. ferus and C. dromedarius was observed within Tylopoda. A markedly higher sequence identity between C. dromedarius and the other vertebrate species was observed at the MSTN locus compared to other genes, thus confirming it as a highly conserved target across mammals.


Introduction
Myostatin (MSTN), a TGF-β family member, is a negative regulator of skeletal muscle development and growth in mammals.The gene sequence has been shown to be highly conserved across vertebrate species (McPherron and Lee, 1997).Loss-of-function mutations associated with increased skeletal-muscle mass and other pleiotropic effects have been detected in several species, including cattle (Grobet et al., 1997), mice (Szabò et al., 1998), humans (Schuelke et al., 2004), sheep (Clop et al., 2006), dogs (Mosher et al., 2007), chickens (Ye et al., 2007;McFarland et al., 2007;Yang et al., 2003) and putatively pigs (Stinckens et al., 2008).These evidences support the hypothesis that also MSTN function could be highly conserved among vertebrates.A total of eight replacement changes at the myostatin gene have been reported in humans (Saunders et al., 2006) while in cattle, one of the most extensively surveyed species, 20 different mutations have been found, 14 of them being located in the coding sequence and six in the intronic regions (for a review, see Stinckens et al., 2011;Rodgers and Garikipati, 2008).In humans, the myostatin gene has been reported to be under the effect of strong positive diversifying selection, as supported by the excess of non-synonymous polymorphisms, to an extent previously observed only for a few other genes, like G6PD and the major histocompatibility complex.Moreover, the frequency patterns observed at two amino acid variants suggested the detected selection signature as a recent phenomenon (Saunders et al., 2006).A recent phenomenon has been postulated also for cattle, where homozygous segments associated to the different MSTN alleles resulted to be longer than those observed in the wild-type counterpart; indeed, the estimate that these alleles had a recent common ancestor less than 400 years ago, together with the evidence of their increasing frequency or fixation have been interpreted as a proof of recent anthropic selection in cattle.A selection pressure on MSTN has been postulated also in bovid where, instead, the observed accelerated evolution of the gene sequence has been explained as the consequence of remote naturally occurring selection (Tellgren et al., 2004).Finally, putative naturally occurring selection has been evoked in salmonids (Ostbye et al., 2007).Based on the above, we decided to partially sequence the MSTN gene in an economical relevant livestock species, Camelus dromedarius, for which only a 256 bp region in the first exon of the MSTN locus was available up to now.Moreover, a preliminary sequence polymorphism analysis was carried out in a population sample of 22 animals from three different Northern African geographic regions.

Sample collection
A total of 22 Camelus dromedarius blood samples from Tunisia (5) Egypt (7) and Algeria (10) were collected from different farms in different areas in order to increase the probability to capture a higher genetic variability at the MSTN locus.

DNA extraction from blood and muscle
Genomic DNA was extracted from blood with FlexiGene DNA Kit (Qiagen), according to the manufacturer's instructions.Briefly 2 ml of blood were mixed with buffer FG1 and centrifuged for 5 min at 2000g.After supernatant was discarded and buffer FG2 /protease mixture was added.Samples were incubated at 65°C for 10 min.After incubation, isopropanol was added for DNA precipitation and samples were centrifuged for 3 min at 2000g.Then ethanol was added and samples were centrifuged for 3 min at 2000g.Supernatant was eliminated and DNA was dissolved with buffer FG3 and incubated for 1 h at 65°C in a water bath.Genomic DNA concentration was evaluated with the Nanodrop 2000C spectrophotometer (Thermo Scientific).Samples were stored a -20°C until use.

PCR amplification
Heterologous oligonucleotide primers to amplify the Camelus dromedarius MSTN locus were designed as follows.First, BLAST analysis using the Vicugna pacos MSTN sequence, available at the Ensemble database (www.ensembl.org)was carried out, which returned, for the Camelus ferus species, a contig (GenBank Accession No AGVR01040332.1)and a scaffold (GenBank Accession No NW_006211358.1),both including the whole-gene MSTN sequence.Alignment of the above sequences was carried out using ClustalW (www.genome.jp/tools/clustalw/) in order to in silico detect possible variant sites.Primers were then designed on the Camelus ferus consensus sequence using Primer3 (www.primer3.ut.ee/, version 4.0.0).Primers and PCR conditions for each target region are provided below.
Exon 1 -The sequence of the forward primer was 5'CCTTGGCATTACTCAAAAGCAA3'; the sequence of the reverse primer was 5'CTTACATACACACGAACAGCC3'. The PCR thermal profile consisted of an initial step at 95°C for 15 min, followed by 35 cycles of 30 s at 94°C, 90 s at 57°C, 90 s at 72°C; the final extension step was carried out at 72°C for 10 min.The sequence also covers putatively 58 bp of the 5'UTR region and 91 bp of the first intron.
Intron 1 (partial) -The sequence of the forward primer was 5'TTTGGGCATTTGCTGACAACCTAGAATGA C3'; the sequence of the reverse primer was 5'CTTTGGGAATCTGAGTAGTTACACTTACT G3'.The PCR thermal profile consisted of an initial step at 95°C for 15 min, followed by 35 cycles of 30 s at 94°C, 90 s at 61°C, 90 s at 72°C; the final extension step was carried out at 72°C for 10 min.
Exon 2 -The sequence of the forward primer was 5'ATCAGATAGTCCTGGAGTAGATCTGCCTT A3'; the sequence of the reverse primer was 5'GTCTGGCTTATGAGCATAGGAAACTGGTA G3'.The PCR thermal profile consisted of an initial step at 95°C for 15 min, followed by 35 cycles of 30 s at 94°C, 90 s at 61°C, 90 s at 72°C; the final extension step was carried out at 72°C for 10 min.The sequence also covers putatively 717 bp of the first intron and 111 bp of the second intron.
Intron 2 (partial) -The sequence of the forward primer was 5'TTACTCTTCTGGCTTACCTGTGCTTGTGT3' ; the sequence of the reverse primer was 5'AAATTTGTCACATCACCAACTTTAATGAT T3'.The PCR thermal profile consisted of an initial step at 95°C for 15 min, followed by 35 cycles of 30 s at 94°C, 90 s at 61°C, 90 s at 72°C; the final extension step was carried out at 72°C for 10 min.
Exon 3 -The sequence of the forward primer was 5'AGCCATAAGGGAAGAATCAAGCC3'; the sequence of the reverse primer was 5'GATCTCCTAGCTCGAAAACTTAGG3'.The PCR thermal profile consisted of an initial step at 95°C for 15 min, followed by 35 cycles of 30 s at 94°C, 90 s at 58°C, 90 s at 72°C; the final extension step was carried out at 72°C for 10 min.The sequence also covers putatively 154 bp of the second intron and 212 bp of the 3'-UTR region.
For all the samples, DNA amplification was carried out in a reaction mixture at a final volume of 25 μl, consisting of Multiplex PCR Master Mix 1X (Qiagen), oligonucleotide primers 0.2 μM (MWG), DNA template 2 ng/μl, and RNAse freewater.

Amplicon purification and sequence analysis
PCR products were separated on a 2% agarose gel electrophoresis at 90V for 1 hour.The fragments were then excised from the agarose gel and purified by using the QIAquick Gel Extraction Kit (Qiagen) according to the manufacturer's instructions.Briefly, DNA was excised from the agarose gel and 3 volumes of buffer QG were added to 1 volume of gel.Gel slices were incubated at 50°C for 10 min.After incubation, isopropanol was added to the sample and the mix was transferred to the column and centrifuged for 1 min at 10000g.For each column 0.5 ml of buffer QG was added and centrifuged for 1 min at 10000g.A washing step was carried out adding 0.75 ml of buffer PE and centrifuging for 1 min at 10000g.The flow-through was discarded and the column was centrifuged for an additional 1 min at 10000g.DNA elution was carried out adding 30 μl of water to the membrane and centrifuging the column 1 min at 10000g.The concentration of the purified PCR products was evaluated with the Nanodrop 2000C spectrophotometer (Thermo Scientific).Sanger sequencing of the purified fragments was performed by the Eurofins MWG Operon service.
Detection of sequence polymorphisms and differentially fixed nucleotide sites was carried out through alignment of Camelus dromedarius sequences against the Camelus ferus contig sequence (GenBank Accession No AGVR01040332.1).Both sequence alignment and pair-wise sequence identity percentages among Camelus dromedarius and other species belonging to the Cetartiodactyla clade were obtained, for the MSTN locus, using ClustalW2 (Larkin et al., 2007).The software Splitstree (Huson and Bryant, 2006) was used to construct a NJ tree from the pair-wise sequence identity matrix.

Results and Discussion
In the present study we report for the first time the nucleotide sequence of the three exons of the MSTN locus for the Camelus dromedarius species obtained after in vitro amplification of genomic DNA using heterologous oligonucleotide primers designed from the publicly available sequence of the Vicugna pacos myostatin gene.Previously, only a 256 bp region in the first exon of the Camelus dromedarius MSTN locus had been published (Shah et al., 2006).Here we sequenced more than 3.6 kb of nucleotide sequence, including the three exons, part of intron 1 and intron 2 and part of the 3' and 5' ends of the Camelus dromedarius myostatin gene.The NJ tree obtained from the pairwise sequence identity among Camelus dromedarius and other species belonging to the Cetartiodactyla clade is provided in Figure 1.Since coding and non-coding regions of a genome evolve at different rates and alignments of non-coding regions are often difficult to interpret, we used here pair-wise sequence identity values estimated after removal of non-exonic sequences.The tree topology was found to be sound and in line with previous reports such that of a closer relationship between Tylopoda and Suiformes previously described, among others, by Kacskovics et al. (2006).As what concerns Tylopoda, a close similarity between C. ferus and C. dromedarius was observed, with a sequence identity of 99.73%,only lower to that observed among Orcinus orca and Tursiops truncatu (99.91%).As expected, a slightly lower identity than that observed between C. dromedarius and C. ferus (Old World camelids) was observed between C. dromedarius and the New World camelid, Vicugna pacos (Table 2).The lowest value was observed between C. dromedarius vs. Bubalus bubalis (93.59%), followed by Bos taurus (93.68%),Capra hircus and Ovis aries (94.04% both), thus supporting the evidence that, despite being generally assimilated to ruminants, camelids have peculiar evolutionary (and physiological) features.Interestingly, when compared to the available literature (Dahiya et al., 2014;Nagarajan et al., 2012Nagarajan et al., , 2011;;Du et al., 2009;Kacskovics et al., 2006), the sequence identity values observed for the MSTN gene between C. dromedarius and the other vertebrate species looks markedly higher, thus confirming previous evidences of the myostatin gene as a highly conserved gene across mammals.Further, to characterize the Camelus dromedarius MSTN gene in terms of sequence polymorphism, we sequenced a total of 22 animals from three different geographic regions (Algeria, Tunisia and Egypt).The results of the polymorphism analysis are represented in Table 1.As can be observed, we only detected three variant nucleotide sites, notably two transitions (798_G/A and 799_C/T) and one transversion (486_G/C).All the polymorphic sites were located in the first intron.Shah et al. (2006) screened a 256 bp region in the first exon of the Camelus dromedarius MSTN locus in 12 samples from six different Pakistani breeds without observing any sequence polymorphism.On the contrary, they found a C→T transition responsible for an Alanine to Valine amino acidic substitution when a 422 bp region of the first exon of the MYF5 gene was screened in the same animals.Vaccarelli et al. (2012) detected, at the T-cell receptor (TCR) gamma SNPs a frequency of one every 114 to 485 nucleotides, depending on the considered rearrangements, and somatic hypermutation was suggested on that locus as a mechanism to diversify the Camelus dromedarius TCR gene repertoire face to the wide spectrum of antigenic determinants presented by a variety of diverse and evolving pathogens.High nucleotide sequence diversity as a consequence of somatic hypermutation had been previously reported also for the T-cell receptor delta chains (Antonacci et al., 2011).We did not observe any polymorphism in the second exon of the myostatin, which is known from other species to harbor functional mutations (Baron et al., 2012).The low diversity observed at the MSTN locus in Camelus dromedarius may reflect the evolutionary history of this species, which likely developed as domesticates from a low variable wild ancestor population.Indeed, a recent world-wide survey of genetic variability in the Camelus dromedarius species carried out adopting both maternally inherited mtDNA and bi-parentally inherited STR markers (Pamela Burger, unpublished data) revealed, consistently with previous published data from STR markers, low levels of genetic variation compared to other livestock species, which has been supposed to derive from the limited geographical distribution of the wild ancestor on the Arabian Peninsula and to the brief co-existence of wild and early domesticated individuals.The extensive dispersal of dromedary camels out of the Arabian Peninsula and the recurrent gene flow phenomena at ancient trading centers along the routes connecting the Arabian Peninsula with the Sahara would also explain the substantial homogeneity observed among our Northern African samples.

Conclusions
This work represented the first extensive characterization of the myostatin gene in the dromedary camel.Interestingly, a very low diversity was observed at the MSTN locus in our Northern African population sample, which may reflect the peculiar evolutionary history of this species.Further ongoing studies from recent available whole genome sequences on the three Camelus species (C.bactrianus, C. ferus and C. dromedarius) will help to better clarify the MSTN patterns of evolution within the interesting Tylopoda genus.

Figure 1 .
Figure 1.NJ tree obtained from pair-wise sequence identity at the three myostatin exons among Camelus dromedarius and other species belonging to the Cetartiodactyla clade.

Table 1 .
Results from the polymorphism analysis at the MSTN locus carried out on 22 Camelus dromedarius samples from three different geographic regions (Algeria, Tunisia and Egypt).Positions are defined in relation to the first nucleotide in the forward primer.

Table 2 .
Matrix of pair-wise percentage sequence identity estimated for the three Camelus dromedarius MSTN exons.