Adaptive molecular evolution of MC1R gene reveals the evidence for positive diversifying selection in indigenous goat populations

Abstract Detecting signatures of selection can provide a new insight into the mechanism of contemporary breeding and artificial selection and further reveal the causal genes associated to the phenotypic variation. However, the signatures of selection on genes entailing for profitable traits between Chinese commercial and indigenous goats have been poorly interpreted. We noticed footprints of positive selection at MC1R gene containing SNPs genotyped in five Chinese native goat breeds. An experimental distribution of F ST was built based on approximations of F ST for each SNP across five breeds. We identified selection using the high F ST outlier method and found that MC1R candidate gene show evidence of positive selection. Furthermore, adaptive selection pressure on specific codons was determined using different codon based on maximum‐likelihood methods; signature of positive selection in mammalian MC1R was explored in individual codons. Evolutionary analyses were inferred under maximum likelihood models, the HyPhy package implemented in the DATAMONKEY Web Server. The results of codon selection displayed positive diversifying selection at the sites were mainly involved in development of genetic variations in coat color in various mammalian species. Positive diversifying selection inferred with recent evolutionary changes in domesticated goat MC1R provides new insights that the gene evolution may have been modulated by domestication events in goats.

transmembrane G-protein-coupled receptor, expressed on melanocyte surface (Javanmard, Arafnajad, Arpanahi, & Moradi, 2015), and it is affecting the production regulation of two pigments; eumelanin and pheomelanin (Liu et al., 2016;Oetting, Austin, & Bennett, 2009). The function of MC1R determines the pigmentation as increased MC1R activity enhances the eumelanin (black/dark, brown pigment) production, whereas decreased MC1R activity enhances pheomelanin (red/ yellow pigment) production (Barsh, 1996;Robbins et al., 1993). MC1R is a single exon gene included 954 bp complete coding sequences with a part of 5′-and 3′-untranslated regions and has seven transmembrane domains with standard 317 amino acid length (Fontanesi et al., 2009;Shimada, Sato, Aplin, & Suzuki, 2009). MC1R gene mutations associated with color phenotype have extensively been studied in different livestock species including goat (Fontanesi et al., 2009), sheep (Våge et al., 2003), and cattle (Klungland, Vage, Gomez-Raya, Adalsteinsson, & Lien, 1995). In goats, a wide variation in coat color has been accounted by agouti loci in several breeds (Robbins et al., 1993). The persistence of dominant black allele and a recessive red allele has been documented in few breeds (Sponenberg, 1990), and on the other hand, wild-type allele of agouti made different phenotypic effects in other species, as in Boer goat the missense mutation (p.k226E) in MC1R gene was associated with red color phenotype (Wu et al., 2006). The MC1R has primary role in coat color evolution; it may be possible to elucidate the long-term coat color evolutionary propensity in a specific ancestry by evaluating the amino acid substitution rate in coding region across the taxa over evolutionary time (Nadeau, Burke, & Mundy, 2007). The identification of genes undergone positive selection has vital role to understand the livestock adaptation to various environmental changes to reveal insights on the selection history. The adaptive evolution was measured by synonymous and nonsynonymous substitution (dN/dS) ratio, as in lion tamarin and mouse, the MC1R sequences have evolved with amino acid changes, which have higher ratio of dN/dS reflecting specific evolutionary episode related coat color (Mundy & Kelly, 2003). Acquisition of new functions in genes is credited to adaptive selection pressures in close association with phenotypes and fitness of organisms (Wright & Rausher, 2010).
Adaptive selective pressure on genes has also been reported to be indicators of functional adaptations developed during the evolution of species that has the tendency of promoting species functional diversification (Vamathevan et al., 2008). Color variations within species present an opportunity to understand the genetic system underlying phenotypic changes and evolutionary meaning of phenotypic characteristics (Hoekstra, Drumm, & Nachman, 2004), as well as coat color provides precious evidence to the natural history of the decree of taxonomic issues for subspecies discrimination (Nachman, Hoekstra, & D'Agostino, 2003). Understanding the evolutionary footprint of MC1R gene will therefore provide valuable information for reconstructing evolutionary history of species and other functional genes, and may provide useful insights into the design of marker-assisted selection and breeding for genetic improvement in goats. In this study, we investigated the molecular evolutionary signatures that may exert selection processes in the MC1R gene in goats and identified evolution footprints that may influence adaptation to different environments.

| Ethics statement
All the essential experimental protocols included in this study were validated by the Law of Animal Husbandry in People's Republic of China (Dec 29, 2005). The entire protocols for sample collection were reviewed and legalized by the Biological Studies, Animal Care and Use Committee of National Animal Husbandry Services, Hubei, China. All measures were adopted to diminish discomfort of the animals during sample collection.

| Animals selection
A total of 526 individuals present in different regions, including Macheng Hubei black-boned goat from Jinyang farm (n = 167), Tonshan Hubei black-boned goat from Tongshan black-boned goat corporation (n = 107), Youzhou black goat from Youzhou black goat farm (n = 138), Nanjing yellow goat (n = 84), and Hubei Boer goat from householders (n = 30) were selected for this study. All goats included in this study were reared under the traditional system and were checked in order to identify domestic variants according to coat color.
The selected animals were unrelated and were from different geographical regions to maximize the genetic variation among individuals.
The sampling was carried out in such a way that maximized spatially unrelated individuals to ensure the breed representativeness and to maximize the genetic diversity among individuals. The genomic DNA samples of 526 individuals were pooled together for selection analysis.

| Sample collection and DNA extraction
Genomic DNA was extracted from blood samples using the TIAGEN DNA Mini Kit (TIAGEN Bio, China). A Nanodrop 2000 (Thermo Scientific, USA) was used to quantify DNA concentration and then stored at −20°C for genotyping.

| Genotyping and polymorphism
The genomic DNA was used to study 13 gene extensions to find out the possible loci in all genomic samples. PCR was carried out on a PCR thermo cycler (Bio-Rad Life Science research) in a total volume of 20 μl containing 12.5 μl PCR master mix (Takara Clontech, Japan), 1.6 μl of both forward and reverse primers, 0.5 μl DNA template, and 5.4 μl nuclease free water. The PCR profile was 35 cycles of 30 s at 95°C, 40 s at 65°C, 40 s at 72°C, 10 min at 72°C and was held at 4°C. The PCR product was analyzed using 2% agarose gel electrophoresis. The PCR products were sequenced using ABI 3100 automated sequencer (Applied Biosystem, Foster city, CA, USA). The primers used in this study are provided in Table S2 (Capra et al., 2017).

| Positive selection analysis
To detect loci under selection, Beaumont and Nichols's approach (Narum & Hess, 2011) implemented in LOSITAN was used. Fixation index (F ST ) and p-values for each locus were calculated using allele frequencies based on heterozygosity simulations, and predictable F ST value was 0.102 in 526 individuals of five populations, and 13 loci. This method presented the contrary selection by finding the outlier with F ST values more than probable overcoming for heterozygosity (Huber, DeGiorgio, Hellmann, & Nielsen, 2016). To build up the population data set, 100,000 simulations were utilized in real data. Quantiles were assumed for transitional joint allocation of F ST against mean heterozygosity with a 95% confidence limit. The outlier was predicted by the loci expressing outside the simulated neutral distribution with differentiating behavior (Antao, Lopes, Lopes, Beja-Pereira, & Luikart, 2008) method, and internal fixed effect likelihood (IFEL) methods Pond & Muse, 2005b;Pond et al., 2006). p-values <.05 were considered as significant. routine, which ensures the robustness against model misspecification over predefined sites through approximate Bayesian method (Murrell et al., 2012;Pond et al., 2011). The fitting of MEME to alignment, MG94xREV codon model, was applied using parameter estimates ω = β/α fitted to the data using the GTR nucleotide model as initial values. The selective pressure was measured with two parameters β: β − ≤ α and β + Yang, 2000), and the alternative model includes four parameters for each site: β − , β + , q − and α estimating site to site substitution variability rates (Pond & Muse, 2005a). The values of p < .05 were considered as significant from the LRT based on χ 2 asymptotic distribution .

| RESULTS
The sequence picture among five goat populations of various breeds revealed thirteen single nucleotide polymorphic sites with an average density 566 bases for each. These single nucleotide polymorphic sites were the consequence of outlier determination. Fixation selection Index-based methods recognized MC1R loci as under selection prance in investigating goat breeds.

| Positive selection of MC1R gene by FDIST analysis
The high F ST outlier method corresponding to F ST distribution was  (Table 2).

| Evolutionary analysis of diversifying selection
The mixed-effect model evolution analysis detected six sites undergone episodic diversifying selection. Among these sites, 59, 86, 108, 153, 179, and 364 were detected under episodic diversifying having p values<.05 (Table 3). This model also estimates the synonymous (α) and nonsynonymous (β) substitution rates, and the sites having values β > α were considered as significant and determined these sites under diversifying selection. The sites 59 was inferred to have experienced pervasive nonsynonymous substitution throughout the evolutionary history with p values .015, and this site evolved with β + > α, as it is under positive selection in all analyses providing the strong evidence of positive selection of MC1R gene, whereas all other sites were conserved. The value of q was determined using Simes' method to reduce the false discovery rate under strict neutral null model (Table 3) (Table 4).

| Evolutionary finger printing of MC1R gene
Evolutionary finger printing of MC1R gene was inferred through codon model evolution based on synonymous and nonsynonymous substitution rates using genetic algorithm. The codon model evolution Pond, Scheffler, Gravenor, Poon, & Frost, 2010) used phylogenetic Markov model, including substitution rates, character frequencies (Pond, Delport, Muse, & Scheffler, 2010), clustering of amino acid substitution rates (Stanfel, 1996), and branch lengths  (Table 5), and these three class rates were further analyzed by genetic algorithm multirate models to evaluate the amino acid substitution rates during evolution. The genetic algorithm models identified the exchange rates at each class rate by estimating the nonsynonymous substitution rates residues labeled by Stanfel class (Figure 3; Stanfel, 1996). The evolutionary rate clusters, indicating the   ing the neutral evolution and only a few sites under the circle above diagonal has positive evolution at the MC1R gene ( Figure 4).

| DISCUSSION
The current advances in the vital index of genetic differences have projected narrative suggestions for exploration of the positive selection goals, which eventually would be supportive to enlighten the genetic drifts and selection roles in evolutionary mechanisms.
Moreover, positive selection signature impedes the genomic regions that play significant roles. Consequently, exploring such regions will provide considerable support for identification of genetic deviations, which would facilitate the distraction of these functional regions and progression in phenotypic assortments. The illuminations of the genetic bases of different traits in most species have been studied by candidate gene approach. The identification of these candidate genes plays an important role in phenotypic variation in livestock population and provides new advancement in the evolutionary process and positive selection (Brown et al., 2013;Storz, Payseur, & Nachman, 2004).
The study of melanin pigmentation in vertebrates has led to the reporting of important target genes that may cause phenotypic variations and divergences in natural populations (Ducrest, Keller, & Roulin, 2008;Hoekstra & Coyne, 2007). One persistent result evolving from most studies is the contribution of the MC1R coding region in explanation of melanism variation, sometimes displaying shared mutations due to convergent evolution between coldly related species (Rompler et al., 2006). MC1R has been shown to play a significant role in various mechanisms such as sexual selection (Mundy et al., 2004), crypsis (Mullen & Hoekstra, 2008), and possibly immunity (Gangoso et al., 2011) although it is generally considered to have few pleiotropic effects (Ducrest et al., 2008). However, in our study, we found molecular evolution in MC1R coding sequences, as would be predictable as they were associated with positively selected regions in nearby locations. Our results are instead consistent with those obtained by Cheviron, Hackett, & Brumfield, 2006). In addition to MC1R, Agouti gene also interacts with MC3R and MC4R and has pleiotropic effects on energy expenditure, food intake (Ducrest et al., 2008).  (2006) showed that using Beaumont and Nichols (Beaumont & Nichols, 1996) approach, it is possible to find candidate gene under selection.
In an integrative approach to detect positive selection, all methods mapped codon evolution at various sites through SLAC, FEL, IFEL, and random effect likelihood and random effect likelihood methods in HyPhy package (Pond & Muse, 2005b). REL detected five codons at sites 59, 218, 248, 327, and 364, while IFEL detected three codons at sites 59, 279, and 327 under positive selection (Table 2). Among these sites, 59 and 327 were detected as positively selected sites in both F I G U R E 4 Evolutionary fingerprints of MC1R gene inferred from alignments on log scale, with the diagonal line corresponding the values α = β for neutral evolution. Dots in the circle indicating the ratio β/α and the area of circle represent the weight of rate classes. The points above the diagonal correspond positive selection and below the diagonal negative selection of the method and the sites 218, 248, and 364 were detected only in REL, it is because other approaches may lack the power of data set containing sequences with low divergence . respectively (Li et al., 2011), as in another study, SLAC detected no site under selection, while FEL and REL detected two sites under positive selection and showed less false-positive results (Sorhannus & Pond, 2006  . Our results contribute the discussion for Bayes factor-based analysis and counting-based method to detect positive selection. REL considered as preferable models to detect the codon sites under positive selection or negative selection (Suzuki, 2004;Suzuki & Nei, 2001).
The mixed-effect model evolution approach based on empirical Bayes procedure identified diversifying selection at MC1R gene, as it estimated the LRT at a particular branch site and identified various sites under diversifying selection. MEME allows site to site and branch to branch distribution of ω and analyze the adaptive evolution at gene level (Yang & Dos Reis, 2011). In our study, we identified six codon sites in the MC1R gene under diversifying selection, and among these sites, there were two sites 59 and 364 (Table 3), which were identified under the selection in former analyses (REL, FEL, SLAC, and IFEL), which were strong evidence of adaptive evolution and selection at MC1R gene. In this analysis of diversifying selection, we restricted the analysis for only those sites with p-value <.05 and empirical Bayes factor threshold > 20, and the best results are achieved when selected branches are placed in conserved lineage (Murrell et al., 2012 (Chen & Sun, 2011;Murrell et al., 2012). Inferring the distribution of genespecific selection parameters potentially improved to detect selection across many sites, and we further performed the FUBAR analysis to manifest the positive selection at MC1R gene. This study identified three sites (Table 4)  The synonymous and nonsynonymous substitutions were identified by codon model substitution over nucleotide and amino acid substitutions (Miyazawa, 2011), and the rate variation significance over protein sites has been validated in nucleotide and amino acid substitution models (Yang, 1994 with maximum-likelihood estimation, which are the computed evolutionary evidence ratio for the gene (Pond, Mannino, Gravenor, Muse, & Frost, 2007). We further performed the evolutionary fingerprint analysis based on probability distribution of site to site synonymous (α) and nonsynonymous (β) substitution rates in alignment and exploited the positive evolution of MC1R gene by estimating the ratio ω = β/α to manifest positive selection, as this ratio (ω = β/α) is an indicator of positive or negative selection in large adaptive evolutionary studies (Arbiza, Dopazo, & Dopazo, 2006;Nielsen et al., 2005).
The MC1R selection patterns discussed here reveal the utility of candidate gene approaches to adaptive phenotypic evolution that can be compared with methods based on genomic scans using neutral markers. Association studies using MC1R and other pigmentation loci may be valuable when studying both Mendelian and quantitative phenotypic traits, either on their own or in coincidence with neutral markers (Nachman, 2005). The drawbacks of candidate gene strategies are obvious when negative results are obtained; that is, there is no association between phenotype and genotype at the candidate locus. The negative inferences in such cases only narrate to the part of the gene studied; for example, in the Phylloscopus study stated by MacDougall-Shackleton et al. (MacDougall-Shackleton, Blanchard, & Gibbs, 2003), it is possible that there is an association with the small portion of MC1R coding region that was not sequenced or with the MC1R promoter region (Mundy, 2005). In general, the success of the approaches depends on various factors such as the accessibility of candidate genes; the evaluation of genetic variation viability in them and the number of loci involved in a particular trait evolution. Fortified selection of candidate genes rendering to the specificity of their phenotypic effects seems to be important, and deleterious effects at loci recognized in inbred captive lines may not always be expressed in captivity.

| CONCLUSIONS
Collection of results from many selection analyses at MC1R gene provides an ideal opportunity to investigate adaptive selection that has influenced the coat color genetic variability and architecture of mammalian genome. As the strong selection of MC1R gene provides the comprehensive map of coat color genetic variability as well as pattern of selection at MC1R gene within the local goat breeds and its correlation among other animal species to reveal molecular evolution. Furthermore, this study will help to construct the data for understanding the selection signature for positive selection of gene among populations in animal breeding and allows the researchers to perform systematic approaches to detect the evolutionary footprints of selection at specific gene.

ACKNOWLEDGMENTS
The author is thankful to anonymous, for their valuable comments and suggestions and proof reading of the manuscript.