# Refining the Camelus dromedarius Myostatin Gene Polymorphism through Worldwide Whole-Genome Sequencing

Silvia Bruno, Vincenzo Landi, Gabriele Senczuk, Samantha Ann Brooks, Faisal Almathen, Bernard Faye, Suheil Semir Bechir Gaouar, Mohammed Piro, Kwan Suk Kim, Xavier David, et al.

<figure><img src="image_1.png" /></figure>

## To cite this version:

Silvia Bruno, Vincenzo Landi, Gabriele Senczuk, Samantha Ann Brooks, Faisal Almathen, et al. Refining the Camelus dromedarius Myostatin Gene Polymorphism through Worldwide Whole-Genome Sequencing. Animals, 2022, 12 (16), pp.2068. 10.3390/ani12162068. [https://doi.org/10.3390/ani12162068](https://doi.org/10.3390/ani12162068)

HAL Id: hal-03772746

https://hal.inrae.fr/hal-03772746

Submitted on 8 Sep 2022

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_2.png" /></figure>
<figure><img src="image_3.png" /></figure>

Article

# Refining the *Camelus dromedarius* Myostatin Gene Polymorphism through Worldwide Whole-Genome Sequencing

Silvia Bruno<fnref n="1" />, Vincenzo Landi<fnref n="2" /> ID, Gabriele Senczuk<fnref n="3" /> ID, Samantha Ann Brooks<fnref n="4" /> ID, Faisal Almathen<fnref n="5" />,<fnref n="6" />, Bernard Faye<fnref n="7" /> ID, Suheil Semir Bechir Gaouar<fnref n="8" />, Mohammed Piro<fnref n="9" />, Kwan Suk Kim<fnref n="10" /> ID, Xavier David<fnref n="11" />, André Eggen<fnref n="11" />, Pamela Burger<fnref n="12" /> ID and Elena Ciani <fnref n="1" />,* ID

<figure><img src="image_4.png" /></figure>

check for
updates

**Citation:** Bruno, S.; Landi, V.; Senczuk, G.; Brooks, S.A.; Almathen, F.; Faye, B.; Gaouar, S.S.B.; Piro, M.; Kim, K.S.; David, X.; et al. Refining the *Camelus dromedarius* Myostatin Gene Polymorphism through Worldwide Whole-Genome Sequencing. *Animals* **2022**, *12*, 2068. https://doi.org/10.3390/ani12162068

Academic Editor: Changxi Li

Received: 30 June 2022

Accepted: 11 August 2022

Published: 14 August 2022

**Publisher's Note:** MDPI stays neutral with regard to jurisdictional claims in published maps and institutional affiliations.

<figure><img src="image_5.png" /></figure>

**Copyright:** © 2022 by the authors. Licensee MDPI, Basel, Switzerland. This article is an open access article distributed under the terms and conditions of the Creative Commons Attribution (CC BY) license (https://creativecommons.org/licenses/by/4.0/).

**Simple Summary:** The dromedary is a multipurpose livestock species. Growing interest in dromedary meat production, increasing intensification of meat production systems and rising business turnover for the dromedary athletic competition sector point to the need for a more effective genetic selection. To this aim, the identification of genetic markers, to be implemented in phenotype–genotype association studies, represents a pre-requisite. The myostatin gene is known to affect, to a various extent, muscularity in several animal species. In this study, we investigated the dromedary myostatin gene sequence variation through the analysis of data from 183 worldwide animals. A total of 99 variants were detected in the target region. Through a bioinformatic approach, we explored the possible functional effects of the detected variants. Several hints emerged, suggesting that natural variants at the dromedary myostatin locus may determine quali-quantitative changes in the myostatin expression, thus likely impacting on dromedary muscularity. Further efforts are needed to collect reliable phenotypic data on dromedary muscularity and racing performances, thus allowing for the unveiling of those genetic markers that are associated with the phenotypic variability in muscularity-related traits. Under intensive rearing systems, large-scale screening of genotypes at those markers in dromedary populations will effectively orient selection decisions and positively impact on rates of genetic gain.

**Abstract:** Myostatin (MSTN) is a highly conserved negative regulator of skeletal muscle in mammals. Inactivating mutations results in a hyper-muscularity phenotype known as “double muscling” in several livestock and model species. In *Camelus dromedarius*, the gene structure organization and the sequence polymorphisms have been previously investigated, using Sanger and Next-Generation Sequencing technologies on a limited number of animals. Here, we carried out a follow-up study with the aim to further expand our knowledge about the sequence polymorphisms at the myostatin locus, through the whole-genome sequencing data of 183 samples representative of the geographical distribution range for this species. We focused our polymorphism analysis on the ±5 kb upstream and downstream region of the MSTN gene. A total of 99 variants (77 Single Nucleotide Polymorphisms and 22 indels) were observed. These were mainly located in intergenic and intronic regions, with only six synonymous Single Nucleotide Polymorphisms in exons. A sequence comparative analysis among the three species within the *Camelus* genus confirmed the expected higher genetic distance of *C. dromedarius* from the wild and

<footnote n="1">¹ Department of Biosciences, Biotechnologies and Biopharmaceutics, University of Bari “Aldo Moro”, 70126 Bari, Italy</footnote>

<footnote n="2">² Department of Veterinary Medicine, University of Bari “Aldo Moro”, Valenzano, 70010 Bari, Italy</footnote>

<footnote n="3">³ Department of Agricultural, Environmental and Food Sciences, University of Molise, 86100 Campobasso, Italy</footnote>

<footnote n="4">⁴ Department of Animal Sciences, University of Florida, Gainesville, FL 32610, USA</footnote>

<footnote n="5">⁵ Department of Public Health, College of Veterinary Medicine, King Faisal University, Al-Ahsa 31982, Saudi Arabia</footnote>

<footnote n="6">⁶ Camel Research Center, King Faisal University, Al-Ahsa 31982, Saudi Arabia</footnote>

<footnote n="7">⁷ CIRAD-ES, UMR SELMET, 34398 Montpellier, France</footnote>

<footnote n="8">⁸ Department of Biology, Abou Bakr Belkaid University of Tlemcen, Tlemcen 13000, Algeria</footnote>

<footnote n="9">⁹ Department of Medicine, Surgery and Reproduction, Institut Agronomique et Vétérinaire Hassan II, Rabat BP 6202, Morocco</footnote>

<footnote n="10">¹⁰ Department of Animal Sciences, Chungbuk National University, Chungbuk 28644, Korea</footnote>

<footnote n="11">¹¹ Illumina, Agrigenomics, 91000 Evry, France</footnote>

<footnote n="12">¹² Research Institute of Wildlife Ecology, Vetmeduni, 1160 Vienna, Austria</footnote>

<footnote>* Correspondence: elena.ciani@uniba.it</footnote>
domestic two-humped camels compared to the genetic distance between *C. bactrianus* and *C. ferus*. In silico functional prediction highlighted: (i) 213 differential putative transcription factor-binding sites, out of which 41 relative to transcription factors, with known literature evidence supporting their involvement in muscle metabolism and/or muscle development; and (ii) a number of variants potentially disrupting the canonical *MSTN* splicing elements, out of which two are discussed here for their potential ability to generate a prematurely truncated (inactive) form of the protein. The distribution of the considered variants in the studied cohort is discussed in light of the peculiar evolutionary history of this species and the hypothesis that extremely high muscularity, associated with a homozygous condition for mutated (inactivating) alleles at the myostatin locus, may represent, in arid desert conditions, a clear metabolic disadvantage, emphasizing the thermoregulatory and water availability challenges typical of these habitats.

**Keywords:** myostatin; dromedary; single nucleotide polymorphisms; indels

# 1. Introduction

Myostatin (*MSTN*), also known as Growth Differentiation Factor-8 (*GDF8*), is a member of the transforming growth factor (TGF-β) superfamily. It is a negative regulator of skeletal muscle development and homeostasis in mammals [1,2]. The gene sequence has been highly conserved across vertebrate species throughout evolution [3]. It consists of three exons and two introns. Myostatin is synthesized as a 375 amino-acid precursor protein called pre-myostatin, that undergoes post-translational modifications in order to become biologically active [4].

Inactivation, through C-terminal deletion, of the *MSTN* gene was shown to determine a two- to three-fold increase in the skeletal muscle mass in mutant compared to wild-type mice [5]. The *MSTN*-knocked out rabbits were shown to generate a heritable phenotype with hyperplasia and hypertrophy of the muscle fibers, which could be valuable in improving rabbits' meat production [6]. MSTN is, indeed, among the most studied genes in meat-producing animals. The mutations at this gene have been shown to naturally occur in several livestock species. In particular, loss-of-function mutations have been associated with a hyper-muscular phenotype, known in cattle as "double muscling" (DBM). The Belgian Blue cattle breed represents a remarkable example of *MSTN*-dependent muscular hypertrophy. In this breed, an 11-bp deletion (g.821-831del11) results in the loss of three amino-acids and a frame-shift mutation, leading to a premature stop codon [7]. The double-muscled cattle display not only an increase in skeletal muscle mass but also a reduction in fat mass [8], with the consequent production of very lean meat. This phenomenon is not surprising, given the fact that *MSTN* has been demonstrated to play key roles not only in myogenesis, but also in adipogenesis [8,9]. In a quantitative trait locus (QTL) survey for growth traits, performed using an F2 Duroc x Pietrain pig resource population, with the Belgian Pietrain breed known for being heavily muscled, Choi et al. [10] identified a locus associated with the *longissimus muscle* area in a region that includes *MSTN*. In the Belgian Texel (Beltex) sheep, Clop et al. [11] reported a mutation creating a potentially illegitimate microRNA target site in the myostatin gene affecting muscularity in sheep. All of the above examples nicely highlight how the Belgian selection of superior-performing animals for extreme muscularity and lean meat percentage has impacted the sequence variability at the *MSTN* locus in these species. Although the double-muscled phenotypes have never been described in the equine species, in thoroughbred horses a SINE insertion in the promoter of the *MSTN* gene inducing a decrease in myostatin expression [12,13] was found to be associated with increased muscle mass [14], thus highlighting that regulatory mutations may more finely tune the hyper-muscularity in livestock species. In addition, the horse SINE also positively affected the aerobic capacity, speed and stamina [15]. A similar phenomenon of increased athletic performance was observed in racing whippet dogs presenting a MSTN mutation leading to a premature stop codon [16].
The dromedary camel (*Camelus dromedarius*) represents an economically relevant livestock species, especially in the countries affected by severe climate change and desertification, assuring the production of milk, meat, leather, wool and a workforce, even in very harsh conditions. The interest in dromedaries and their products is also rapidly growing in Europe where camel farms have been established, such as “La Camélerie” (https://www.lacamelerie.fr/ accessed on 7 August 2022), in France, and the Kamelenmelkerij Smits (https://www.kamelenmelk.nl/ accessed on 7 August 2022), in Netherlands. Moreover, the Emirates Industry for Camel Milk and Products (EICMP), owner of the “Camelicious” brand (https://camelicious.com/ accessed on 7 August 2022), has become the first Middle Eastern company to obtain a license to export dromedary products to the European Union. In recent years, there has been an increase in the demand for camel meat as an alternative lean food, characterized by a lower fat content, lower cholesterol and a higher proportion of polyunsaturated fatty acids than most consumed meat from cattle and sheep [17]. In addition, the dromedary is traditionally exploited, mainly in the Arabian Peninsula countries, for racing and beauty contests. Based on the above, its morphological features of muscularity may influence the economic value of the animal. Nevertheless, despite the fact that differences in muscularity were observed among distinct populations [18] and/or individuals [19], to our knowledge no evident hyper-muscular phenotype has so far been described in dromedaries.

In *C. dromedarius*, the *MSTN* gene-structure organization and the sequence polymorphisms have been previously investigated, using Sanger [20] and Next-Generation Sequencing (NGS) [21] technologies. No report of loss-of-function mutations exists so far in this species, and the above mentioned studies highlighted a lower degree of sequence variation than in other livestock species, possibly related to its peculiar evolutionary history, characterized by: (i) a domestication occurring from a markedly bottlenecked wild genetic stock [22] and (ii) a potentially not particularly strong human-based positive selection pressure on the muscularity traits. Indeed, although camels were largely used as pack animals in antiquity, the obvious attitudinal advantage deriving from higher muscle mass has been attained preferentially through inter-specific hybridization with *C. bactrianus* [23]. Moreover, the adaptive phenotypic and genetic plasticity phenomena typical of mammals living in xeric environments may well have represented negative constraints on muscle mass development in order to cope with thermoregulatory and water availability challenges [24,25]. In order to check whether the low sequence variation previously reported [20,21] in *C. dromedarius* at the *MSTN* locus could also be affected by the relatively limited number of animals enrolled in the studies, we extended the NGS analysis here to a world-wide set of 183 dromedaries. The sequence polymorphisms results are presented, and their possible functional effect is predicted through in silico analysis.

**2. Materials and Methods**

**2.1. *Sampling and Whole-Genome Sequencing***

Within the frame of an Illumina®-funded project (2019 Agricultural Greater Good Initiative), through a large international collaborative effort, 161 *Camelus dromedarius* biological samples (whole blood and hair follicles) were collected from 18 countries (Table S1, Supplementary Materials), representative of the geographic distribution range of this species. Whenever possible, the samples were collected to avoid closely unrelated animals, based on herders' knowledge. The vast majority of the samples were collected for previous research projects (FWF P24706-B25 and P29623-B25, PI: P. Burger; EU Arimnet2 CARAVAN, PI: E. Ciani) during routine veterinary procedures, therefore no ethical assessment and permit was specifically required for sampling (art. 1, comma 5 of the Directive 2010-63-EU). Notwithstanding that, procedures to avoid any possible animal distress during the samples' collection were put in place, according to local consolidated practices. The hair samples from Ethiopia were specifically collected for this project by the researchers of Ethiopian Biotechnology Institute, with no ethical authorization required for hair sampling (art. 1, comma 5 of the Directive 2010-63-EU). In order
to assure local cultures' sensitivity [26] in the traditional dromedary countries, communication tools based on dialogue as a basis for trust-building and informed (oral) consent were adopted.

The transfer of biological material was completed following the example material transfer agreement provided as the annex to the FAO guideline on the molecular genetic characterization of animal genetic resources (FAO, 2011).

The DNA from whole blood was extracted using the Qiagen DNeasy Blood and Tissue kit (Qiagen, Venlo, The Netherlands). The DNA from the hair follicles was extracted with a DNA salting-out method [27]. The DNA quantity and quality was assessed with a Qubit fluorometer (Invitrogen, Vienna, Austria).

Paired-end library preparation and sequencing steps were carried out by the Illumina Solutions Center in San Diego, CA. In brief, the DNA was first quantified using the Quant-iT PicoGreen dsDNA Assay Kit (Thermo Fisher Scientific, Waltham, MA, USA). Between 100–500 ng of DNA was used as the input to prepare the libraries, using the Illumina DNA Prep library kit and the IDT for Illumina DNA/RNA UD Index set (Illumina, San Diego, CA, USA), following the reference guide. The final libraries were quantified using the Quant-iT PicoGreen dsDNA Assay Kit and normalized to 2 nM using an average region size of 600 bps. The libraries with unique barcodes were pooled together and sequenced on the Illumina NovaSeq 6000 system using the NovaSeq 6000 S4 300 Cycle Reagent Kits (Illumina, San Diego, CA, USA). In addition to the above, we also capitalized the publicly available whole-genome sequence data for this species (Table S2, Supplementary Materials), by retrieving 22 sequence sets representative of nine countries.

## 2.2. *Mapping and Variant Calling*

The raw reads from the 183 considered animals were analyzed by using the DRAGEN Germline App v.3.9.5 on the BaseSpace™ platform (Illumina, San Diego, CA, USA), using the default setting parameters. The DRAGEN workflow includes both alignment and variant calling algorithms: the reads were first mapped and aligned to the *C. dromedarius* CamDro3 reference genome (Genebank accession: GCA_000803125.3), then the variant calling was performed producing a genotype vcf.gz output file overall, including 17,679,716 variants.

## 2.3. *Characterization of Variants in the Myostatin Locus*

The genotype file was loaded into the RStudio software v.1.4.1103 (R version 4.0.3) and the variants in the range of ±5 kb upstream and downstream of the myostatin gene were retained by using the Tidyverse package v.1.3.2. Overall, our target region was 16,757 bps long, including 6757 bps of the *MSTN* gene and 5000 bps upstream and downstream (chr5: 58,454,553–58,471,309). The variants were classified based on their nature (SNPs or indels), their localization (5'-UTR, exonic, intronic, 3'-UTR, intergenic) and their occurrence in the previously published paper by Favia et al. [21]. In addition, the popular [28] Tajima's D neutrality test [29], as implemented in VCFtools v0.1.13 [30], was carried out on the chromosome 5 harboring the *MSTN* locus, performing the analysis on 2000 bps windows. If the D values are too large or too small, the neutral null hypothesis is rejected. As a rule of thumb, we considered here that values greater than +2 (balancing selection) or less than -2 (directional selection) are likely to be significant [28]. For the SNPs detected in the exonic regions, the effect of the polymorphisms on the *MSTN* amino-acid sequence and the codon usage frequency patterns were evaluated. The latter were assessed via the ATGme software [31], using the *C. dromedarius* codon usage table available at http://www.kazusa.or.jp/codon/cgi-bin/showcodon.cgi?species=9838 (accessed on 7 August 2022).

## 2.4. *Inter-Specific Comparative Sequence Variant Analysis*

A sequence comparative analysis of the target myostatin region described above was carried out by BLASTing, the publicly available *Camelus dromedarius* myostatin reverse-complement sequence (GenBank accession: GCA_000803125.3) against *Camelus bactrianus* (GenBank accession: GCA_000767855.1) and *Camelus ferus* (GenBank accession: GCA_009834535.1) correspond-
ing sequences. We opted here for using the reverse-complement sequence of the *C. dromedarius* CamDro3 reference genome (Genebank accession: GCA_000803125.3) in order to allow for an easier comparison with the previously published paper by Favia et al. [21]. Hence, all of the subsequent analyses were also performed on the CamDro3 reverse-complement sequence.

## 2.5. *In Silico Functional Prediction*

The web-based software TFBIND [32], available at https://tfbind.hgc.jp/ (accessed on 7 August 2022), was used to identify the transcription factor binding sites (TFBS) and their possible disruption due to the presence of Single Nucleotide Polymorphism or indels. The TFBIND analysis was carried out, considering the 2 kb upstream the *C. dromedarius* myostatin gene. In addition, the online bioinformatic tool Human Splicing Finder (HSF) [33], available at https://www.genomnis.com/ (accessed on 7 August 2022), was used to perform a Sequence Analysis (Version 1.5.1.). The HSF tool employs several matrices to predict the variants' putative effect on the splicing motifs, including the acceptor and donor splice sites, the branch point and auxiliary sequences known to either enhance or repress the splicing. Both of the above analyses were carried out, for each variant, using the two input sequences harboring the alternative alleles.

# 3. Results and Discussion

## 3.1. Sequence Analysis

Our sequence variant analysis performed in the ±5 kb upstream and downstream region of the *MSTN* gene highlighted the presence of a total of 99 variants (Table S3, Supplementary Materials).

Out of them, 77 were Single Nucleotide Polymorphisms (SNPs) and 22 were indels. Based on the localization, 72 variants were intergenic (namely, 38 variants in the 5 kb region upstream of the *MSTN* gene and 34 variants in the 5 kb region downstream of the *MSTN* gene), 1 SNP in the 5' UTR region, 1 SNP in exon 1, 11 variants in intron 1, 1 SNP in exon 2, 9 variants in intron 2, and 4 SNPs in exon 3 (Table S3, Supplementary Materials). The 11 SNPs, previously detected by Favia et al. [21], in the sequence overlapping with our target region using a dataset of nine dromedaries, were also found in our 183 samples' dataset. Out of the six SNPs that we detected in the exonic regions of the myostatin gene, none of them was responsible for the amino acid substitution at the protein level (Table 1).

<table><caption>**Table 1.** Analysis of potential amino acid substitutions caused by SNPs in the coding region.</caption><thead><tr><th rowspan="2">#<br>Exon</th><th rowspan="2">Variant<br>Name</th><th rowspan="2">Project<br>ID *</th><th rowspan="2">Chromosome<br>Position **</th><th colspan="2">Allele</th><th rowspan="2">Base Codon<br>Position</th><th rowspan="2">Mutation<br>Type</th><th rowspan="2">Allele<br>Frequency</th></tr><tr><th>Allele 1</th><th>Allele 2</th></tr></thead><tbody><tr><td>1</td><td>5var58465820</td><td>PRJEB55295</td><td>58465820</td><td>T</td><td>C</td><td>3rd</td><td>Silent</td><td>0.5</td></tr><tr><td>2</td><td>5var58463788</td><td>PRJEB55295</td><td>58463788</td><td>G</td><td>T</td><td>3rd</td><td>Silent</td><td>0.5</td></tr><tr><td>3</td><td>5var58461341</td><td>PRJEB55295</td><td>58461341</td><td>G</td><td>A</td><td>3rd</td><td>Silent</td><td>0.5</td></tr><tr><td>3</td><td>5var58461338</td><td>PRJEB55295</td><td>58461338</td><td>A</td><td>T</td><td>3rd</td><td>Silent</td><td>0.5</td></tr><tr><td>3</td><td>5var58461332</td><td>PRJEB55295</td><td>58461332</td><td>T</td><td>C</td><td>3rd</td><td>Silent</td><td>0.5</td></tr><tr><td>3</td><td>5var58461263</td><td>PRJEB55295</td><td>58461263</td><td>A</td><td>G</td><td>3rd</td><td>Silent</td><td>0.5</td></tr></tbody></table>

The frequency of codon usage was tested, for the six SNPs in Table 1. Two rare codons (CCG and GTA, Allele 2 codons in Table S4, Supplementary Materials) were detected in the reference *MSTN* sequence in exon 1 (frequency 8.1‰) and exon 2 (frequency 5.5‰). These were converted into common codons (CCA and GTC, Allele 1 codons in Table S4, Supplementary Materials) when the alternative alleles detected in our study were fed into the ATGme software. On the other hand, when considering the alternative alleles, a codon in exon 3 (CAT) resulted in being rare (6.4‰), while being common when considering the reference sequence (codon CAC).

An average density of one SNP every 217 bps was observed in our 16,757 bps considered region. This value is much higher than that reported by Favia et al. [21], who mention

<footnote>* Accession number for the project submitted to European Variant Archive (EVA, https://www.ebi.ac.uk/eva/ accessed on 7 August 2022). ** The position refers to the CamDro3 assembly (GCA_000803125.3).</footnote>
an average density of one SNP every 1.5 kbps. It must be noted that the number of whole-genome sequenced animals differed markedly between the two studies, with nine animals sampled from the Gulf countries (Kingdom of Saudi Arabia, United Arab Emirates, Qatar), Pakistan, Sudan, Kenya and Spain (Canary islands) in Favia et al. [21], and 183 animals representative of 19 countries from three continents, in this study, respectively. Despite the fact that the average SNP density observed in this study is generally comparable with the expected figures in mammals (e.g., in the ovine species, on average one SNP every 204 bps, by re-sequencing the nine animals drawn from different breeds [34]); in the bovine species, on average one SNP every 104 or 434 bps by sequencing the leptin or the amyloid precursor protein genes, respectively, in 22 individuals from the two subspecies *Bos taurus* and *Bos indicus* [35]) and in aves (e.g., in the chicken, about one SNP every 200 bps, by re-sequencing three animals from a corresponding number of domestic chicken breeds contrasted with the sequence of their wild ancestor, the red jungle fowl [36]); in ducks, on average one SNP about every 86 bps, through the genotyping-by-sequencing of 49 animals from the same flock [37]), still, the fact that, to reach such a value, a number of animals much higher than in the above mentioned studies had to be sequenced points out that a generally lower sequence variability exists in *C. dromedarius*, thus confirming our previous results on a more limited numbers of animals [20,21].

A rather clear conservative evolutionary constraint appears to exist on the myostatin locus, as highlighted by the much higher level of polymorphism in the intergenic region (73%) compared to that in the transcribed region, and the lack of missense/nonsense/frame-shift mutations in the exons. This phenomenon is not uncommon. It has been reported, for the human species, that, on average, about 69% of the SNPs are located in intergenic regions versus 31% of the SNPs located in genic regions [38]. Out of them only a small fraction (0.76%) has been reported to be located in the exons and, roughly, half of them have been described as missense mutations, while a 3:1 proportion of the missense to silent mutations would be expected under neutrality (absence of selection) [38]. In the *C. dromedarius* *MSTN* gene, the proportion of exonic versus genic SNPs was 6:27 (roughly, 0.22%), with no missense or nonsense SNPs observed, possibly reflecting an evolutionary constraint. Indeed, under the harsh desert conditions, more than in other environments, hyper-muscularity may represent a clear metabolic disadvantage, emphasizing thermoregulatory and water availability challenges typical of these habitats. A support to this hypothesis is also offered by the allele frequency patterns observed for the six exonic SNPs detected at the *C. dromedarius* *MSTN* gene, where the balanced allele frequencies (Table 1) nicely match with the results of the Tajima's D neutrality test (Table S5, Supplementary Materials), the latter pointing to a possible balancing selection effect.

## 3.2. *Comparative Analysis*

A comparative analysis carried out by BLASTing the *C. dromedarius* myostatin sequence (GCA_000803125.3) against the corresponding sequences in *C. ferus* (GCA_009834535.1) and *C. bactrianus* (GCA_000767855.1), for the 99 loci segregating within the *C. dromedarius* species in the region, including the ±5 kbp upstream and downstream of *MSTN*, highlighted that *C. bactrianus* and *C. ferus* displayed the same allele (i.e., the same nucleotide) of the *C. dromedarius* reference sequence at 36 and 48 loci, respectively. Out of them, 31 and 43 loci, respectively, harbored the allele shown to be the most frequent within our *C. dromedarius* considered cohort (consensus allele, in Table S6, Supplementary Materials). The nucleotide at the G/T SNP locus (position 58,466,681 bps) segregating within *C. dromedarius* displays, in *C. ferus*, a single base deletion. The number of loci where *C. bactrianus* and *C. ferus* share the same allele is 85/99, in line with the results from previous studies on divergence time evaluation within the *Camelus* genus [39,40], thus pointing, as expected, to a genetic distance of *C. dromedarius* from the wild and domestic two-humped camels higher than the genetic distance between *C. bactrianus* and *C. ferus*.
### 3.3. *In Silico* Functional Prediction

A bioinformatic analysis using the TFBIND tool was conducted to identify whether the allelic variations of the SNPs and indels in the target region resulted in differential binding of any transcription factors (TFs). Overall, in the 2 kbps proximal to the transcription initiation site of the *C. dromedarius* myostatin gene, 11 variants (i.e., nine SNPs and two indels) were identified. By using the two alternative alleles for each variant, a total of 213 differential putative TF-binding sites (TFBS) (out of which 82 were unique) were observed (Table S7, Supplementary Materials); out of which 41 (14 being unique) relative, overall, to ten TFs with known literature evidence supporting their involvement in muscle metabolism and/or muscle development, i.e., COMP1 (cooperates with myogenic proteins 1, 2), C/EBP (CCAAT/enhancer binding protein, 12), MYOG (Myogenin, 2), COUP (Chicken Ovalbumin Upstream Promoter-transcription factor II, 2), MEF2 (Myocyte Enhancer Factor-2, 16), SRF (Serum Response Factor, 1), DELTAEF1 (Zinc Finger E-Box Binding Homeobox 1, 1), PBX1 (PBX Homeobox 1, 3), MYOD (Myoblast Determination Protein 1, 1) and E2F (Retinoblastoma-Binding Protein, 1). Out of them, MEF2, DELTAEF1 and MYOG are involved in promoting and regulating the skeletal muscle cell differentiation program during myogenesis [41–43]. Moreover, MYOG and MYOD belong to a family of proteins known as myogenic regulatory factor (MRF) that plays a major role in regulating the skeletal muscle differentiation [44]. Interestingly, MYOD may also regulate muscle repair [45]. SRF is required for MyoD expression and, consequently, for the myogenic differentiation and maintenance of muscle fibers [46]. A study conducted in *Xenopus laevis* suggested that the RSRFC4 transcription factor may regulate muscle-specific transcription in embryos, and may have other roles during muscle development [47]. Besides, in the skeletal muscles, E2F is essential for the full activation of the myogenic genes during muscle growth and myofibrillogenesis [48], and COMP1 plays a role in the interaction with myogenin [49]. PBX1 is considered a pioneer factor, since it is constitutively bound to the promoter of the myogenin gene, serving as a platform for the MYOD binding in inactive chromatin, thereby preparing the genes of the skeletal muscle lineage for activation [50]. Lastly, COUP-TFII inhibits myogenesis and skeletal muscle metabolism in vitro and in vivo, by repressing the myoblast fusion [51]. Similarly, all of the C/EBP isoforms have a role in the inhibition of myogenic differentiation [52,53].

We next focused our attention on the possible mechanism of alternative splicing; therefore, the HSF tool was used to predict the effects of the identified variant sites on splicing signals. In detail, HSF was designed to perform in silico predictions for the formation or disruption of splice donor sites, splice acceptor sites, branching points (BPs), exonic splicing silencer sites (ESSs) and exonic splicing enhancer sites (ESEs). We performed the mentioned analysis on the variants falling in the exonic, intronic and UTR regions of the gene by separately feeding HSF, for each variant, with the sequences harboring the two alternative alleles. As for the effects on the splice donor and acceptor sites (Table 2), out of the 27 analyzed variants (20 intronic, 6 exonic and 1 in the 5' UTR), six variants were predicted to disrupt three natural splice sites and to create four alternative splice sites (three acceptor and one donor splice sites).

<table><caption>**Table 2.** Differential acceptor and donor splice sites identified by the HSF sequence analysis using the two sequences harboring the alternative alleles.</caption><thead><tr><th colspan="4">ALLELE 1</th><th colspan="3">ALLELE 2</th></tr><tr><th>INTRONIC VARIANTS</th><th>Type <sup>1</sup></th><th>Motif <sup>2</sup></th><th>Value <sup>3</sup></th><th>Type <sup>1</sup></th><th>Motif <sup>2</sup></th><th>Value <sup>3</sup></th></tr></thead><tbody><tr><td>5var58464953</td><td>-</td><td>-</td><td>-</td><td>Acceptor splice site</td><td>CCCGTCTGCAGAT</td><td>81.96</td></tr><tr><td>5var58464910</td><td>Donor splice site</td><td>TATGTTATT</td><td>72.09</td><td>-</td><td>-</td><td>-</td></tr><tr><td>5var58464775</td><td>Acceptor splice site</td><td>GCACCTTAACAGAG</td><td>77.89</td><td>-</td><td>-</td><td>-</td></tr><tr><td>5var58464463</td><td>-</td><td>-</td><td>-</td><td>Donor splice site</td><td>GGTGTTAAT</td><td>66.52</td></tr><tr><td>5var58464142</td><td>Donor splice site</td><td>GAAGTAGGT</td><td>81.61</td><td>Acceptor splice site</td><td>GAAAGGAAGCAGGT</td><td>68.36</td></tr><tr><th>EXONIC VARIANTS</th><th>Type <sup>1</sup></th><th>Motif <sup>2</sup></th><th>Value <sup>3</sup></th><th>Type <sup>1</sup></th><th>Motif <sup>2</sup></th><th>Value <sup>3</sup></th></tr><tr><td>5var58461332</td><td>-</td><td>-</td><td>-</td><td>Acceptor splice site</td><td>ATTGCACCTAAGAG</td><td>74.03</td></tr></tbody></table>

<footnote n="1">¹ Type of splicing signals: acceptor splice sites or donor splice sites; ² Sequence motif harboring the predicted splicing signals; ³ The signals values range from 65 (weak) to 100 (very strong).</footnote>
Concerning the effects on the branching points, eight SNPs were predicted to cause the creation of eight putative branching points and the disruption of three putative BPs (Table S8, Supplementary Materials). A total of 22 variants were found to alter the ESEs (Table S9, Supplementary Materials) and 25 variants were found to alter the ESSs (Table S10, Supplementary Materials). In particular, as to what concerns the variants responsible for the formation/disruption of the splice donor/acceptor sites, we focused our attention on: (i) the single T/TAATAA indel at position 58,464,910 bps, as this was the only indel among those predicted to alter the splice donor/acceptor sites and (ii) the SNP at position 58464775, as this was the only locus to be responsible for the alteration of all of the categories of the splice elements (donor/acceptor sites, branching points, exonic splicing silencer/enhancer sites). Simulating retention of the 793 bps region of intron 1, subsequent to the creation of an alternative donor splice site when the TAATAA allele at position 58,464,910 bps is present, the expected protein would be prematurely truncated due to the presence of a stop codon only two amino acids downstream of the canonical end of the peptide encoded by the exon 1 (Figure 1A). Given the potential relevance of such an observation, we checked whether the TAATAA allele was present in a homozygous state in at least one animal in our cohort, but it appeared that the seven copies of the TAATAA allele were all in a heterozygous state in the seven animals harboring the insertion. These animals were sampled in different countries (Tunisia, Libya, Mauritania, Sudan), with the exception of two samples coming from Morocco but belonging to two different sub-populations (Guerzni and Marmouri). Although all of the animals bearing the TAATAA allele were of African (mainly North-African) origin, still, the geographically scattered distribution of this allele may suggest the presence of homozygous animals as rather uncommon. Similarly to our results, Ceccobelli et al., also in Marchigiana cattle, [54] did not found homozygote-mutated animals in the investigated dataset of 78 samples, and observed only a light superiority in the live weight reached by heterozygous bulls compared to wild-type bulls, although it was not significant. It is indeed known that the incredibly high muscularity typical of the double muscle phenotype generally occurs in association with a homozygous condition for the mutated allele at the myostatin locus. Mosher et al. [16] discovered a 2-bp deletion in the whippet dog *MSTN* gene that, in the homozygote state, results in a double-muscling phenotype commonly referred to as the “bully” whippet, while the heterozygote animals display a phenotype of intermediate musculature. Qian et al. [55] demonstrated that the homozygote *MSTN* mutant pigs had an apparent double muscle phenotype, and individual muscle mass increased by 100% over their wild-type controls at eight months of age as a result of myofiber hyperplasia, while the *MSTN*-mutant heterozygotes animals had much lower muscle mass than the homozygote *MSTN* mutant pigs.

Interestingly, a prematurely truncated protein would also be expected when simulating the retention of the 871-bps region of intron 1, subsequent to the creation of an alternative acceptor splice site when the C allele at position 58,464,775 bps is present instead of the G allele (Figure 1B). A total of 78 animals harbored, in our study cohort, the C allele and, intriguingly, all of them presented a heterozygous genotype. This time, almost all of the countries represented in our 183 sample set had at least one heterozygous animal, with Morocco having the highest number of heterozygous animals (13).

The alternative splicing mechanisms can affect gene expression, thereby generating several distinct proteins [56] and enriching the phenotypic diversity among the individuals in a population [57]. To date, few studies focusing on the on-lab detection of alternative splicing have been conducted in *C. dromedarius*, among which are included the one by Premraj et al. [58] that reported four tissue-specific transcript variants of the interleukin-26 (IL-26) gene; the one by Ryskaliyeva et al. [59] that identified two isoforms of the dromedary αs2-CN milk phosphoprotein/protein arising from differential splicing events, increasing the ability of camel caseins to generate potentially bioactive peptides; and the one by Kappeler et al. [60] that highlighted the presence of two lactophorin variants, resulting from alternative mRNA splicing, in lactating mammary glands of dromedaries. As for the
myostatin gene, a previous study from our group [21] did not find any evidence for the
alternative splicing in this species.

(A)

MQLQIYVYIYLFMLIVAGPVDLNENNEQKENVEKEGLCNACMWRQNTKSSRLEAIKIQILSKLRLE
TAPNISKDAIRQLLPKAPPLRELIDQYDVQRDDSSDGSLEDDDYHATTETIITMPTECKSY
CISTILLTVVLVFMRNRSIFRLFQAVRVYVSRREKSFIFQDFMRNLMRLRAAALYIPGSKAKNKK
NACIALILYSLAPICLCFHSLIPPRQGLDTANNVKNTNLQQLVITNMHLQQTFILITEKGRSL
MLKPIWAFADNLELLLFKTISHSVFMFFTNYNIFKLLHWKVKYFKIICISKSNDEVNIAKS
ANDKMPRNKHFNQIGYGSQSPTYTLAMVLLLRIPGLQISRQAHCLITFN

Methionine (start codon) in red; exon 1 highlighted in yellow, stop codons highlighted in red.

(B)

MQLQIYVYIYLFMLIVAGPVDLNENNEQKENVEKEGLCNACMWRQNTKSSRLEAIKIQILSKLRLE
TAPNISKDAIRQLLPKAPPLRELIDQYDVQRDDSSDGSLEDDDYHATTETIITMPTEKA
QERFLCHVSVILLYP LTGHKTVSKYCLSYVIRLCNKFSFKTVLSDSPGVDLPYLIMENQKERNF
VFLLTMTACP RQCSLTFEIA MRHLNFGIN LLVLFLSKAIPAPKHNRCTIFSTNSR
LLSVLSPGNSGLGEGFLLPD
LVQLLSKCNYSDSORILSGCSSTVSLVLSNHHHFKISSTVHSFIEFSLVHSSLERKQVPHKQLKKNRA
KNSKQLLFCLIFFVMNEILHSFSLKKTNKYLLFNRKMLTDYGGVCSFSKIISITLFFSYSFI

Methionine (start codon) in red; exon 1 highlighted in yellow, stop codons highlighted in red.

**Figure 1.** Expected amino acid sequences based on HSF results. (A) Simulation of retention of 793 bps in intron 1 due to the creation of an alternative donor splice site when the TAATAA allele belonging to 5var58464910 is present; (B) Simulation of retention of 871 bps in intron 1 due to the creation of an alternative acceptor splice site when the C allele belonging to 5var58464775 is present.

**4. Conclusions**

This study contributed to obtaining a wider and clearer picture of the sequence varia-
tion at the *C. dromedarius* myostatin locus. Previously undescribed polymorphisms (SNPs
and indels) were detected, and the variant sites previously published by our group were
confirmed in a larger world-wide cohort. The myostatin gene was shown to harbor a low
sequence variability in *C. dromedarius*, possibly reflecting a specific evolutionary constraint.

On the other hand, the environmental and/or human-mediated selection toward inter-
mediately high muscular mass, that is more compatible (i.e., less metabolically challenging)
with life in extreme xeric habitats compared to extreme, homozygous hyper-muscularity,
may be the reason for the observed lack, in our study cohort, of animals homozygous
for the alleles in silico predicted to significantly alter the myostatin-splicing process and,
consequently, its amino-acidic sequence and function.

Further efforts need to be devoted to collecting reliable phenotypic data on dromedary
muscularity and racing performances, thus allowing the unveiling of possible genotype-
phenotype associations useful (i) to better understand the molecular bases of muscle-mass
development in dromedaries and (ii) to possibly assist selection for meat and athletic traits
in this species via the large-scale genotyping of validated genetic markers. Indeed, although
in a non-uniform manner, the intensification of dromedary meat production systems is
rapidly ongoing in many of the traditional dromedary countries, together with a growing
business turnover for the dromedary athletic competition sector. Under the above scenario,
it is very likely that genotype-assisted selection may soon also come under the spotlight in
this species.

**Supplementary Materials:** The following supporting information can be downloaded at: https://www.mdpi.com/article/10.3390/ani12162068/s1, Table S1: Details about the 161 *C. dromedarius* animals specifically sampled for this study; Table S2: Details about the 22 publicly available *C. dromedarius* whole-genome sequence data used in this study; Table S3: List of the 99 sequence variants detected in the range ±5 kbps upstream and downstream of the myostatin gene (chr5: 58,454,553–58,471,309); Table S4: Results of the codon usage frequency analysis carried out on the exonic SNPs of the myostatin gene; Table S5: Results of the Tajima's D test carried out on chromosome 5 harboring the *MSTN* gene; Table S6: Results of the inter-specific comparative sequence variant analysis within the *Camelus* genus; Table S7: Results of the TFBS analysis for the sequence variants falling in the 2 kb upstream the *C. dromedarius* myostatin gene; Table S8: Differential HSF Branch Point matrix results obtained by using the two sequences harboring the alternative alleles; Table S9: Differential Exonic Splicing Enhancer (ESE) matrices results obtained by using the two sequences harboring the alternative alleles; Table S10: Differential Exonic Splicing Silencer (ESS) matrices results obtained by using the two sequences harboring the alternative alleles.

**Author Contributions:** Conceptualization, E.C. and P.B.; methodology, E.C., S.B., P.B.; software, X.D. and A.E.; validation, S.B., V.L. and G.S.; formal analysis, S.B., V.L. and G.S.; investigation, E.C. and S.B.; resources, S.A.B., F.A. and B.F.; data curation, S.B.; writing—original draft preparation, S.B. and E.C.; writing—review and editing, E.C., S.B., P.B., X.D., A.E., V.L., G.S., K.S.K., S.A.B., F.A., B.F., S.S.B.G. and M.P.; visualization, S.B.; supervision, E.C.; project administration, E.C.; funding acquisition, E.C., P.B., S.A.B. and F.A. All authors have read and agreed to the published version of the manuscript.

**Funding:** This research was funded by Illumina (United States), grant 2019 Agricultural Greater Good Initiative and by the EU PRIMA Consortium, grant CAMEL-SHIELD.

**Institutional Review Board Statement:** Not applicable.

**Informed Consent Statement:** Not applicable.

**Data Availability Statement:** The variant data for this study have been deposited in the European Variation Archive (EVA) at EMBL-EBI under accession number PRJEB55295 (https://www.ebi.ac.uk/eva/?eva-study=PRJEB55295).

**Conflicts of Interest:** The authors have read the journal's guidelines and have the following competing interests: André Eggen and Xavier David work within the Illumina® company, who provided data, used in this work, generated as part of the 2019 Illumina® Greater Good Initiative. This program spurs critically needed research that will increase the sustainability, productivity and nutritional density of an agriculturally important crop and livestock species. Grant recipients receive donations of Illumina® products to support their projects. The other authors have no competing interests.

**References**

1. Beyer, T.A.; Narimatsu, M.; Weiss, A.; David, L.; Wrana, J.L. The TGFβ Superfamily in Stem Cell Biology and Early Mammalian Embryonic Development. *Biochim. Biophys. Acta* **2013**, *1830*, 2268-2279. [CrossRef] [PubMed]

2. McCoy, J.C.; Walker, R.G.; Murray, N.H.; Thompson, T.B. Crystal Structure of the WFIKKN2 Follistatin Domain Reveals Insight into How It Inhibits Growth Differentiation Factor 8 (GDF8) and GDF11. *J. Biol. Chem.* **2019**, *294*, 6333-6343. [CrossRef] [PubMed]

3. McPherron, A.C.; Lee, S.J. Double Muscling in Cattle Due to Mutations in the Myostatin Gene. *Proc. Natl. Acad. Sci. USA* **1997**, *94*, 12457-12461. [CrossRef] [PubMed]

4. Wolfman, N.M.; McPherron, A.C.; Pappano, W.N.; Davies, M.V.; Song, K.; Tomkinson, K.N.; Wright, J.F.; Zhao, L.; Sebald, S.M.; Greenspan, D.S.; et al. Activation of Latent Myostatin by the BMP-1/Tolloid Family of Metalloproteinases. *Proc. Natl. Acad. Sci. USA* **2003**, *100*, 15842-15846. [CrossRef]

5. McPherron, A.C.; Lawler, A.M.; Lee, S.J. Regulation of Skeletal Muscle Mass in Mice by a New TGF-Beta Superfamily Member. *Nature* **1997**, *387*, 83-90. [CrossRef]

6. Zheng, Y.; Zhang, Y.; Wu, L.; Riaz, H.; Li, Z.; Shi, D.; Rehman, S.; Liu, Q.; Cui, K. Generation of Heritable Prominent Double Muscle Buttock Rabbits via Novel Site Editing of Myostatin Gene Using CRISPR/Cas9 System. *Front. Vet. Sci.* **2022**, *9*, 842074. [CrossRef]

7. Grobet, L.; Martin, L.J.; Poncelet, D.; Pirottin, D.; Brouwers, B.; Riquet, J.; Schoeberlein, A.; Dunner, S.; Ménissier, F.; Massabanda, J.; et al. A Deletion in the Bovine Myostatin Gene Causes the Double-Muscled Phenotype in Cattle. *Nat. Genet.* **1997**, *17*, 71-74. [CrossRef]
8. Deng, B.; Zhang, F.; Wen, J.; Ye, S.; Wang, L.; Yang, Y.; Gong, P.; Jiang, S. The Function of Myostatin in the Regulation of Fat Mass in Mammals. *Nutr. Metab.* **2017**, *14*, 29. [CrossRef]

9. Rebbapragada, A.; Benchabane, H.; Wrana, J.L.; Celeste, A.J.; Attisano, L. Myostatin Signals through a Transforming Growth Factor Beta-like Signaling Pathway to Block Adipogenesis. *Mol. Cell. Biol.* **2003**, *23*, 7230–7242. [CrossRef]

10. Choi, I.; Steibel, J.P.; Bates, R.O.; Raney, N.E.; Rumph, J.M.; Ernst, C.W. Application of Alternative Models to Identify QTL for Growth Traits in an F2 Duroc x Pietrain Pig Resource Population. *BMC Genet.* **2010**, *11*, 97. [CrossRef]

11. Clop, A.; Marcq, F.; Takeda, H.; Pirottin, D.; Tordoir, X.; Bibé, B.; Bouix, J.; Caiment, F.; Elsen, J.-M.; Eychenne, F.; et al. A Mutation Creating a Potential Illegitimate MicroRNA Target Site in the Myostatin Gene Affects Muscularity in Sheep. *Nat. Genet.* **2006**, *38*, 813–818. [CrossRef] [PubMed]

12. Santagostino, M.; Khoriauli, L.; Gamba, R.; Bonuglia, M.; Klipstein, O.; Piras, F.M.; Vella, F.; Russo, A.; Badiale, C.; Mazzagatti, A.; et al. Genome-Wide Evolutionary and Functional Analysis of the Equine Repetitive Element 1: An Insertion in the Myostatin Promoter Affects Gene Expression. *BMC Genet.* **2015**, *16*, 126. [CrossRef]

13. O'Hara, V.; Cowan, A.; Riddell, D.; Massey, C.; Martin, J.; Piercy, R.J. A Highly Prevalent SINE Mutation in the Myostatin (MSTN) Gene Promoter Is Associated with Low Circulating Myostatin Concentration in Thoroughbred Racehorses. *Sci. Rep.* **2021**, *11*, 7916. [CrossRef] [PubMed]

14. Petersen, J.L.; Mickelson, J.R.; Rendahl, A.K.; Valberg, S.J.; Andersson, L.S.; Axelsson, J.; Bailey, E.; Bannasch, D.; Binns, M.M.; Borges, A.S.; et al. Genome-Wide Analysis Reveals Selection for Important Traits in Domestic Horse Breeds. *PLoS Genet.* **2013**, *9*, e1003211. [CrossRef] [PubMed]

15. Rooney, M.F.; Hill, E.W.; Kelly, V.P.; Porter, R.K. The "Speed Gene" Effect of Myostatin Arises in Thoroughbred Horses Due to a Promoter Proximal SINE Insertion. *PLoS ONE* **2018**, *13*, e0205664. [CrossRef]

16. Mosher, D.S.; Quignon, P.; Bustamante, C.D.; Sutter, N.B.; Mellersh, C.S.; Parker, H.G.; Ostrander, E.A. A Mutation in the Myostatin Gene Increases Muscle Mass and Enhances Racing Performance in Heterozygote Dogs. *PLoS Genet.* **2007**, *3*, e79. [CrossRef]

17. Suliman, G.M.; Alowaimer, A.N.; Hussein, E.O.S.; Ali, H.S.; Abdelnour, S.A.; El-Hack, M.E.A.; Swelum, A.A. Chemical Composition and Quality Characteristics of Meat in Three One-Humped Camel (Camelus Dromedarius) Breeds as Affected by Muscle Type and Post-Mortem Storage Period. *Animals* **2019**, *9*, 834. [CrossRef]

18. Abdallah, H.R.; Faye, B. Phenotypic classification of Saudi Arabian camel (*Camelus dromedarius*) by their body measurements. *Emir. J. Food Agric.* **2012**, *24*, 272–280.

19. Kadim, I.T.; Mahgoub, O.; Purchas, R.W. A Review of the Growth, and of the Carcass and Meat Quality Characteristics of the One-Humped Camel (Camelus Dromedaries). *Meat Sci.* **2008**, *80*, 555–569. [CrossRef]

20. Muzzachi, S.; Oulmouden, A.; Cherifi, Y.; Yahyaoui, H.; Zayed, M.; Burger, P.; Lacalandra, G.; Faye, B.; Ciani, E. Sequence and Polymorphism Analysis of the Camel (Camelus Dromedarius) Myostatin Gene. *Emir. J. Food Agric.* **2015**, *27*, 367. [CrossRef]

21. Favia, M.; Fitak, R.; Guerra, L.; Pierri, C.L.; Faye, B.; Oulmouden, A.; Burger, P.A.; Ciani, E. Beyond the Big Five: Investigating Myostatin Structure, Polymorphism and Expression in Camelus Dromedarius. *Front. Genet.* **2019**, *10*, 502. [CrossRef] [PubMed]

22. Almathen, F.; Charruau, P.; Mohandesan, E.; Mwacharo, J.M.; Orozco-terWengel, P.; Pitt, D.; Abdussamad, A.M.; Uerpmann, M.; Uerpmann, H.-P.; De Cupere, B.; et al. Ancient and Modern DNA Reveal Dynamics of Domestication and Cross-Continental Dispersal of the Dromedary. *Proc. Natl. Acad. Sci. USA* **2016**, *113*, 6707–6712. [CrossRef] [PubMed]

23. Berthon, R.; Mashkour, M.; Burger, P.; Çakırlar, C. Domestication, Diffusion and Hybridization of the Bactrian Camel: A Zooarchaeological Perspective. In *Les Vaisseaux du désert et des Steppes: Les Camélidés dans l'Antiquité (Camelus dromedarius et Camelus bactrianus)*; Redon, B., Agut-Labordère, D., Eds.; Archéologie(s); MOM Éditions: Lyon, France, 2020; pp. 1–26.

24. Curtin, N.A.; Bartlam-Brooks, H.L.A.; Hubel, T.Y.; Lowe, J.C.; Gardner-Medwin, A.R.; Bennitt, E.; Amos, S.J.; Lorenc, M.; West, T.G.; Wilson, A.M. Remarkable Muscles, Remarkable Locomotion in Desert-Dwelling Wildebeest. *Nature* **2018**, *563*, 393–396. [CrossRef] [PubMed]

25. Newman, M.T. Biological Adaptation of Man to His Environment: Heat, Cold, Altitude, and Nutrition. *Ann. N. Y. Acad. Sci.* **1961**, *91*, 617–633. [CrossRef]

26. Cooper, T.L.; Kirino, Y.; Alonso, S.; Lindahl, J.; Grace, D. Towards Better-Informed Consent: Research with Livestock-Keepers and Informal Traders in East Africa. *Prev. Vet. Med.* **2016**, *128*, 135–141. [CrossRef] [PubMed]

27. De Volo, S.B.; Reynolds, R.T.; Douglas, M.R.; Antolin, M.F. An Improved Extraction Method to Increase DNA Yield from Molted Feathers. *Condor* **2008**, *110*, 762–766. [CrossRef]

28. Doğan, İ.; Doğan, N. Statistical Tests for Neutrality: Review. *Turk. Klin. J. Biostat.* **2017**, *9*, 167–174. [CrossRef]

29. Tajima, F. Statistical Method for Testing the Neutral Mutation Hypothesis by DNA Polymorphism. *Genetics* **1989**, *123*, 585–595. [CrossRef]

30. Danecek, P.; Auton, A.; Abecasis, G.; Albers, C.A.; Banks, E.; DePristo, M.A.; Handsaker, R.E.; Lunter, G.; Marth, G.T.; Sherry, S.T.; et al. The Variant Call Format and VCFtools. *Bioinformatics* **2011**, *27*, 2156–2158. [CrossRef]

31. Daniel, E.; Onwukwe, G.U.; Wierenga, R.K.; Quaggin, S.E.; Vainio, S.J.; Krause, M. ATGme: Open-Source Web Application for Rare Codon Identification and Custom DNA Sequence Optimization. *BMC Bioinform.* **2015**, *16*, 303. [CrossRef]

32. Tsunoda, T.; Takagi, T. Estimating Transcription Factor Bindability on DNA. *Bioinformatics* **1999**, *15*, 622–630. [CrossRef] [PubMed]

33. Desmet, F.-O.; Hamroun, D.; Lalande, M.; Collod-Béroud, G.; Claustres, M.; Béroud, C. Human Splicing Finder: An Online Bioinformatics Tool to Predict Splicing Signals. *Nucleic Acids Res.* **2009**, *37*, e67. [CrossRef] [PubMed]
34. Kijas, J.W.; Townley, D.; Dalrymple, B.P.; Heaton, M.P.; Maddox, J.F.; McGrath, A.; Wilson, P.; Ingersoll, R.G.; McCulloch, R.; McWilliam, S.; et al. A Genome Wide Survey of SNP Variation Reveals the Genetic Structure of Sheep Breeds. *PLoS ONE* **2009**, *4*, e4668. [CrossRef] [PubMed]

35. Vignal, A.; Milan, D.; SanCristobal, M.; Eggen, A. A Review on SNP and Other Types of Molecular Markers and Their Use in Animal Genetics. *Genet. Sel. Evol.* **2002**, *34*, 275–305. [CrossRef]

36. International Chicken Polymorphism Map Consortium. A Genetic Variation Map for Chicken with 2.8 Million Single-Nucleotide Polymorphisms. *Nature* **2004**, *432*, 717–722. [CrossRef]

37. Zhu, F.; Cui, Q.-Q.; Hou, Z.-C. SNP Discovery and Genotyping Using Genotyping-by-Sequencing in Pekin Ducks. *Sci. Rep.* **2016**, *6*, 36223. [CrossRef]

38. Zhao, Z.; Fu, Y.-X.; Hewett-Emmett, D.; Boerwinkle, E. Investigating Single Nucleotide Polymorphism (SNP) Density in the Human Genome and Its Implications for Molecular Evolution. *Gene* **2003**, *312*, 207–213. [CrossRef]

39. Geraads, D.; Didier, G.; Barr, A.; Reed, D.; Laurin, M. The Fossil Record of Camelids Demonstrates a Late Divergence between Bactrian Camel and Dromedary. *APP* **2020**, *65*, 251–260. [CrossRef]

40. Wu, H.; Guang, X.; Al-Fageeh, M.B.; Cao, J.; Pan, S.; Zhou, H.; Zhang, L.; Abutarboush, M.H.; Xing, Y.; Xie, Z.; et al. Camelid Genomes Reveal Evolution and Adaptation to Desert Environments. *Nat. Commun.* **2014**, *5*, 5188. [CrossRef]

41. Gossett, L.A.; Kelvin, D.J.; Sternberg, E.A.; Olson, E.N. A New Myocyte-Specific Enhancer-Binding Factor That Recognizes a Conserved Element Associated with Multiple Muscle-Specific Genes. *Mol. Cell. Biol.* **1989**, *9*, 5022–5033.

42. Jethanandani, P.; Kramer, R.H. Alpha7 Integrin Expression Is Negatively Regulated by DeltaEF1 during Skeletal Myogenesis. *J. Biol. Chem.* **2005**, *280*, 36037–36046. [CrossRef] [PubMed]

43. Cao, Y.; Kumar, R.M.; Penn, B.H.; Berkes, C.A.; Kooperberg, C.; Boyer, L.A.; Young, R.A.; Tapscott, S.J. Global and Gene-Specific Analyses Show Distinct Roles for Myod and Myog at a Common Set of Promoters. *EMBO J.* **2006**, *25*, 502–511. [CrossRef] [PubMed]

44. Botzenhart, U.U.; Gerlach, R.; Gredes, T.; Rentzsch, I.; Gedrange, T.; Kunert-Keil, C. Expression Rate of Myogenic Regulatory Factors and Muscle Growth Factor after Botulinum Toxin A Injection in the Right Masseter Muscle of Dystrophin Deficient (Mdx) Mice. *Adv. Clin. Exp. Med.* **2019**, *28*, 11–18. [CrossRef] [PubMed]

45. Goudenege, S.; Pisani, D.F.; Wdziekonski, B.; Di Santo, J.P.; Bagnis, C.; Dani, C.; Dechesne, C.A. Enhancement of Myogenic and Muscle Repair Capacities of Human Adipose-Derived Stem Cells with Forced Expression of MyoD. *Mol. Ther.* **2009**, *17*, 1064–1072. [CrossRef] [PubMed]

46. Coletti, D.; Daou, N.; Hassani, M.; Li, Z.; Parlakian, A. Serum Response Factor in Muscle Tissues: From Development to Ageing. *Eur. J. Transl. Myol.* **2016**, *26*, 6008. [CrossRef]

47. Chambers, A.E.; Kotecha, S.; Towers, N.; Mohun, T.J. Muscle-Specific Expression of SRF-Related Genes in the Early Embryo of Xenopus Laevis. *EMBO J.* **1992**, *11*, 4981–4991. [CrossRef]

48. Zappia, M.P.; Rogers, A.; Islam, A.B.M.M.K.; Frolov, M.V. Rbf Activates the Myogenic Transcriptional Program to Promote Skeletal Muscle Differentiation. *Cell Rep.* **2019**, *26*, 702–719.e6. [CrossRef]

49. Funk, W.D.; Wright, W.E. Cyclic amplification and selection of targets for multicomponent complexes: Myogenin interacts with factors recognizing binding sites for basic helix-loop-helix, nuclear factor 1, myocyte-specific enhancer-binding factor 2, and COMP1 factor. *Proc. Natl. Acad. Sci. USA* **1992**, *89*, 9484–9488. [CrossRef]

50. Grebbin, B.M.; Schulte, D. PBX1 as Pioneer Factor: A Case Still Open. *Front. Cell Dev. Biol.* **2017**, *5*, 9. [CrossRef]

51. Lee, H.-J.; Kao, C.-Y.; Lin, S.-C.; Xu, M.; Xie, X.; Tsai, S.Y.; Tsai, M.-J. Dysregulation of Nuclear Receptor COUP-TFII Impairs Skeletal Muscle Development. *Sci. Rep.* **2017**, *7*, 3136. [CrossRef]

52. Marchildon, F.; Lala, N.; Li, G.; St-Louis, C.; Lamothe, D.; Keller, C.; Wiper-Bergeron, N. CCAAT/Enhancer Binding Protein Beta Is Expressed in Satellite Cells and Controls Myogenesis. *Stem Cells* **2012**, *30*, 2619–2630. [CrossRef] [PubMed]

53. Yang, H.; Mammen, J.; Wei, W.; Menconi, M.; Evenson, A.; Fareed, M.; Petkova, V.; Hasselgren, P.-O. Expression and Activity of C/EBPβ and δ Are Upregulated by Dexamethasone in Skeletal Muscle. *J. Cell. Physiol.* **2005**, *204*, 219–226. [CrossRef] [PubMed]

54. Ceccobelli, S.; Perini, F.; Trombetta, M.F.; Tavoletti, S.; Lasagna, E.; Pasquini, M. Effect of Myostatin Gene Mutation on Slaughtering Performance and Meat Quality in Marchigiana Bulls. *Animals* **2022**, *12*, 518. [CrossRef] [PubMed]

55. Qian, L.; Tang, M.; Yang, J.; Wang, Q.; Cai, C.; Jiang, S.; Li, H.; Jiang, K.; Gao, P.; Ma, D.; et al. Targeted Mutations in Myostatin by Zinc-Finger Nucleases Result in Double-Muscled Phenotype in Meishan Pigs. *Sci. Rep.* **2015**, *5*, 14435. [CrossRef] [PubMed]

56. Matlin, A.J.; Clark, F.; Smith, C.W.J. Understanding Alternative Splicing: Towards a Cellular Code. *Nat. Rev. Mol. Cell Biol.* **2005**, *6*, 386–398. [CrossRef]

57. Manning, K.S.; Cooper, T.A. The Roles of RNA Processing in Translating Genotype to Phenotype. *Nat. Rev. Mol. Cell Biol.* **2017**, *18*, 102–114. [CrossRef]

58. Premraj, A.; Nautiyal, B.; Aleyas, A.G.; Rasool, T.J. Identification of Interleukin-26 in the Dromedary Camel (Camelus Dromedar-ius): Evidence of Alternative Splicing and Isolation of Novel Splice Variants. *Mol. Immunol.* **2015**, *67* Pt B, 357–368. [CrossRef]

59. Ryskaliyeva, A.; Henry, C.; Miranda, G.; Faye, B.; Konuspayeva, G.; Martin, P. Alternative Splicing Events Expand Molecular Diversity of Camel CSN1S2 Increasing Its Ability to Generate Potentially Bioactive Peptides. *Sci. Rep.* **2019**, *9*, 5243. [CrossRef]

60. Kappeler, S.; Farah, Z.; Puhan, Z. Alternative Splicing of Lactophorin mRNA from Lactating Mammary Gland of the Camel (Camelus Dromedarius). *J. Dairy Sci.* **1999**, *82*, 2084–2093. [CrossRef]