All sampling procedures in this study was adhered to the guidelines of Ethic Committee of Veterinary Research Institute of Mongolia (ref. N 256/2010). Goats samples were collected across the middle Mongolian region incorporating desert, dry steppe, mountain steppe and dense forest natural zones with decreasing temperature going northward (including partial northern, central and southern Mongolia) (Figure 1). A total of 105 goat samples from three breeds were used for this study; 85 MNG (accounted for >95% of total goats in Mongolia), 10 GGS and 10 ZJT. The GGS and ZJT breeds were introduced from Soviet-bloc countries [19] and have become well adapted in the desert region in southern Mongolia, mostly reared in Bulgan soum (or village) of Omnogobi province, and Jinst soum in Bayankhongor province, respectively. While MNG is the most widely distributed and common across Mongolia. We sampled from eight different populations; six populations for MNG and one unique population for GGS and ZJTs, respectively. The phenotypes (coat color) of these animals are shown in Figure 1. The GGS and ZJT goat breeds’ status were characterized by their uniform black and white coat colors [19]. Whereas, in MNG goats variable coat colors, such as, red and black color, white, blue, white spotted on the body or face or red-headed white, or various combinations of these colors have been observed (Supplementary Table S1).

Eight mL of blood samples were collected from MNG and GGS goats, whereas, tissue samples from ear were collected from ZJT goats. Blood (with vacuum tubes inclduing ethylenediaminetetraacetic acid) and ear (with 70% ethanol) samples were stored in car portable freezer (−4°C to −12°C) during prolonged field work. Blood samples of Mongolian goats were preserved on FTA classic cards (Whatman: Lot No.9767339, Maidstone, UK) in order to transfer samples to Korea. Genomic DNA was extracted from three punches of filter paper (FTA, 2 mm punches) using MagMax DNA Multi-Sample Ultra Kit (TermoFisher Scientific, Baltics UAB V.A. Vilnius, Lithuania) using an extraction protocol with a magnetic beads following the manufacturer’s instructions.

A long exon (938-bp) of the MC1R gene of goat was amplified to identify the polymorphisms. A pair of primers (Forward: 5′-CTG AGA GCA AGC ACC CTT TC -3′; Reverse: 5′-GTC CTG TGA TTC CCC TCT CA -3′) for caprine MC1R polymerase chain reaction (PCR) amplification and sequencing was designed with Primer3Plus tools (online available at: https://primer3plus.com/cgi-bin/dev/primer3plus.cgi) used previously and published as goat reference sequence (Genbank accession number: Y13958) [20]. The PCR amplifications for goat MC1R gene was performed by using thermal cycler C1000 (Bio-Rad laboratories Inc., Hercules, CA, USA) system in a volume of 20 μL that contained 2.0 μL genomic DNA (15 ng/μL), 0.8 μL of each primer (10 pmoles/μL), 2.0 μL 10× reaction buffer, 1.6 μL dNTP, 0.2 μL prime Taq DNA polymerase, and 12.6 μL distilled water for adjusting final PCR master mix volume. The reaction was carried out with initial denaturation for 3 min at 94°C, followed by 35 cycles of denaturation for 30 s at 94°C, annealing for 30 s at 58°C, extension for 45 s at 72°C, and final extension at 72°C for 5 min. The amplified DNA fragments were detected on 1.5% agarose gel, and size checked using 100 bp DNA ladder (GeNetBio, Daejeon, Korea), and photographed under UV light (in Printgraph, Atto Corporation, Tokyo, Japan). For sequencing, PCR amplicons were purified with PCR purification kit (GeNetBio, Korea) following the manufacturer’s recommendation. All 105 samples were sequenced in this study (Cosmogenetech Company [www.cosmogenetech.com]).

A total of 105 novel sequences of Mongolian goats (85 MNG, 10 GGS, and ZJTs, respectively) were initially assembled in BioEdit V7.2.5 [21] and aligned with the Clustal W multiple aligning tools [22] for identification of SNPs and further analysis. The SNPs were identified with the Clustel W aligner and DnaSP v.5.10.01 [23]. The variant effect predictor (VEP) in Ensemble web tool (http://asia.ensembl.org/Tools/VEP) was used to obtain the information of SNPs’ nomenclature and prediction of consequences of variants on the protein sequence (e.g. synonymous, missense and non-sense) based on NCBI accession number Y13958 as reference. The Mongolian goat MC1R genotypes were obtained from analyzing and comparing the obtained sample using a sequencing electropherogram approach (Supplementary Figure S1a).

The DnaSP V5.10.01 (DNA sequence polymorphism, [23]) was applied for estimating the diversity parameters, including haplotype diversities (Hd), and nucleotide diversities (Pi), as well as number of haplotypes (H), number of segregating sites (S), average number of pairwise different (K), and haplotype file were also generated for further haplotype network analysis. The molecular phylogenetic tree was constructed by maximum-likelihood (ML) method based on Tamura-Nei model [24] in molecular evolutionary genetics analyser V7 (MEGA) [25]. We drew the median-Joining (MJ) network using NETWORK V5.0.0.1 to reveal genetic relationships among MC1R haplotypes in Mongolian goats [26].

Alternatively, the online tool PROVEAN V1.1 (Protein Variation Effect Analyzer, avialable at: http://provean.jcvi.org/seq_submit.php) was used to predict amino acid substitution effect on protein function. If the Provean provided a value that was equal to or below a predefined threshold (e.g. −2.5), the protein variant was predicted to have a deleterious effect on its function, while the protein variant was predicted to have a neutral effect when the score was above the threshold [27]. Finally, the PROTTER v1.0 online tools (http://wlab.ethz.ch/protter/#) [28] was used to predict the secondary structure of the caprine MC1R protein using reference amino acid sequence (CAR82366) [3].

The protein sequence for MC1R gene (Genbank accession number: CAR82366) of capra hircus was acquired from National Center for Biotechnology Information (NCBI) and blast in Clustal W tools [22] against protein sequences, that also were obtained from NCBI: ovis aries (SBU12237), bos taurus (NP_776533), equus caballus (NP_001108006), Sus scrofa (ALN40197), Homo sapiens (NP_002377), Felis catus (NP_ 001009324), Canus l. Familiaris (NP_001014304), Oryctolagus cuniculus (CBJ17605), Mus musculus (NP_032585), Gallus gallus (AGY49247), Anas platerhynchos (NP_001297734), Danio rerio (NP_851301).

For the purpose of statistics, coat colors were divided into three main groups: i) Red (n = 50, complete red, white-faced red and red-pied), ii) Black (n = 31, complete black and black-pied), and iii) White (n = 18). These three color groups also are dominant in Mongolian goat populations. Six goats were excluded from association test due to low frequency of other minor coat color expression. In this study, the two main statistical tests were used to ascertain association between genotype and phenotype, a logistic regression for coat color groups of black versus white and a chi-square (x²) test for three groups of coat color, and these tests were implemented in R version 3.5.1.

The PCR amplification produced 1,020-bp fragments of coding region (CDS). After initial sequence editing, the 938-bp long fragments were utilized in downstream analysis. The obtained sequences represented three different breeds from eight independent populations.

Analysis of the MC1R gene sequences revealed a total of five substitutions in Mongolian goat populations. The number of polymorphic sites for MNG, GGS, and ZJTs goats were comparable. In addition, the overall estimated nucleotide (Pi) diversity for investigated goat populations was 9.7×10⁻⁴ (Supplementary Table S2). Alternatively, Pi value was estimated for the three coat color groups. According to diversity estimation, the highest Pi was observed for white animals (11.6×10⁻³) (Supplementary Table S3). Observed diversity indices may suggest that black and white uniformed animals had higher polymorphism than predominant red.

A total of seven MC1R gene haplotypes were defined by five polymorphic sites in Mongolian goats. The similar number of MC1R haplotypes and much higher number of agouti signaling protein haplotypes were found in different goat populations [3,29,30]. The MJ network and ML tree shows their relationships (Figures 2, Supplementary Figure S2). According to these clusters, H2 (TGAGT), H3 (TGATT), and H5 (TGAGA) were closely clustered compared to others (Supplementary Figure S2), while H6 (CGGTT) was separated by three mutation steps (c.748T>G, c.183C>T, and c.676A>G) from central H2 (Figure 2). In addition, Table 1 reports haplotype frequencies and their distributions among Mongolian goat populations. The H2 showed the highest frequency (60.4%) in two breeds of goat (or populations; without ZJTs) including most of the recorded phenotypes (9 coat colors), followed by H3 (21.6%), while relatively low number of animals (from 1 to 6) were associated to the remainder of haplotypes (Table 1, Figure 2). In particularly, H7 (TAATT) acounted for only 1 white animal. Intrestingly, one haplotype (H5) was unique for black goats, whereas the majority of red goats (86.1% of total red group) belonged to the main MC1R haplotypes (H2), white predominant in H3 (55.5%) and H6 (22.2%), and black goats (50.0%) goat predominant in H4 (TAAGT) (Figure 2). According to the findings of this study, unique H5 (TGAGA) may have some relationship with black uniform, unfortunately there were a lack of such animals sample.

The CDS region in goat MC1R gene encodes a protein of 317 amino acids that was similar to bos taurus and ovis aries with 97% and 99% identity, respectively (Figures 3, Supplementary Figure S1b). A total of five SNPs (polymorphic sites) were identified upon aligning the obtained sequences (CDS) of Mongolian goats (Table 2). Two out of five SNPs were identified as synonymous mutations (c.183C>T, c.489G>A), while the remaining three were missense mutations (non-synonymous), c.676A>G, c.748T>G, c.770T>A causing the p.K226E, p.F250V, and p.F257Y amino acid substitutions in transmembrane domains, TM1, TM4, TM6, and intracellular loops-3 of the MC1R protein (Table 2, Figure 3). Meanwhile, lysine (at 226) and phenylalanine (at 250 and 257) were highly conserved in the aligned melocortin 1 receptor in the amino acid sequences of diverse mammalian and avian species (Supplementary Figure S1b).

Among the identified five SNPs, three SNPs were previ ously described, namely, synonymous mutation c.183C>T and missense mutation c.748T>G by Fontanesi et al [3], and c.676A>G by Wu et al [17] and Javanmard et al [31]. Therefore, SNPs c.489G>A and c.770T>A are uniquely observed in this study. The Provean score predicted the amino acids substitution effect on protein function (Table 2). The Provean score revealed that two of three missense mutations (p.F250V and p.F257Y) may have deleterious effects with −6.435 and −2.674 scores, respectively. Whereas, p.K226E mutation had a neutral effect on protein function with −1.566 score (Table 2). Fontanesi et al [3] reported that putative functional effect of missense p.F250V mutation on the protein was deleterious with a subPSEC score −5.998 (substitution position-specific evolutionary conservation). Moreover, the p.F250V seems to be a more deleterious effect, because of this mutation scored a higher functional effect score in this and a previous study [3]. The subPSEC score works in parallel with Provean score to predict functional effect on protein indunction. No functional effect score was avialable for previously described missense p.K226E mutation [17].

A total of 105 animals were genotyped to analyze polymorpism in caprine MC1R gene. In addition, all of these animals’ phenotypes (coat color) were recorded during field sampling. Supplementary Table S4 reveals mutated allele and genotype frequencies for five identified polymorphic sites on the three phenotypic groups in Mongolian goats. Here, we further highlighted two common missense mutations (c.676A>G and c.748T>G), none for their two synonymous, and a very rare missense mutaion (c.770T>A). We performed a chi-sqaure test to assess genotype frequency differences among the three coat color groups, while a logistic regression for association of black versus white phynotypes and genotypes additionally.

For misse nse c.676A>G (p.K226E) mutation, frequencies of genotypes AG and GG was significantly greater in white (x² = 39.1; d.f = 4; p<0.0001), while the AA genotypes are present in both red and black (Supplementary Table S4). In addition, dominant allele A (E^D), frequency also was significantly higher, especially in the red group (see, Supplementary Table S5; Figure S3). In a previous study, Wu et al [17] found a higher frequency of GG genotype for the majority of investigated animals (dominated by complete white), while only Boer (red-headed white) goats carried AA genotype for c.676A>G mutation. Further, the authors highlighted that G allele in this mutation may be a loss of function mutation associated with complete white goat coat colors [17]. Whereas, Javanmard et al [31] found the higher frequency of AA genotype in eumelanin phenotypes.

Whereas, black and red, particularly red group accounted for higher frequencies of mutated allele (G) (84%) than white had for c.748T>G (p.F250V) mutation. However, genotype frequencies varied significantly among the groups, and homozygous genotype GG significantly higher in the red and black (x² = 34.3; d.f = 4; p<0.0001) (Supplementary Table S4). It seems mutated type allele (G, or recessive e) may contribute to the red phenotypes in Mongolian goats (Supplementary Tables S4, S5). In contrast, Fontanesi et al [3] reported that goats with black eyed-white and black or brown coats were found to have a higher frequency of GG genotype (0.94% and 1.0%, respectively), while red goats (Derivata di Siria) had only 0.26% mutated allele for mutation c.748T>G.

Additionally, the main two color variations in mammals that are controlled by mutations in MC1R gene are black and white. A logistic regression also revealed a highly significant values for association between genotypes and these two phenotypes (p<0.05). Detailed information regarding the results between white and black individuals using logistic regression analysis were illustrated in Table 3. Therefore, based on the findings obtained from this study (association tests), the Extension locus in goat might play in some active role in balck versus white coat color varation, at least in the three investigated breeds in this study (white animals from MNG and ZJT; black from MNG and GGS).

Unfortunately, the most rare missense mutation c.770T>A (p.F257Y), which may contribute to the black coat color in Mongolian goats, effects could not be positively determined due to only six complete black goats (animals from both MNG and GGS; mutated allele frequency 0.11%) were heterozygous in this mutation (Supplementary Table S5). Moreover, five of these six goats belong to H5 haplotype (TGAGA). To reach a definitive conclusion, more animals need to be sampled.

Previously, some mutations at the extension locus that af fected different coat colors in different mammalian species such as, in goat [3], sheep [8,31], cattle [9], pigs [10], horses [11], dogs [12], rabbits [14] were reported. Whereas, in goat extension locus, Fontanesi et al [3] reported higly significant associations (Fisher exact test, p<0.001) between amino acid substitution and black coat color and red spotted white goats for c.801C>G missense and c.673T nonsense mutations, respectively. Unfortunately, the aforementioned mutations were not recorded in this study. In particularly, Fontanesi et al [3] emphasized the nonsense mutation (c.673T)’s alteration that resulted for shorter MC1R protein (Figure 3).

Interestingly, in this study two common missense muta tions (c.676A>G and c.748T>G) were expressed in animals with white uniform (did not influence any melanin) tended to carry a higher frequency of heterozygote genotypes (0.72% and 0.66%, respectively), this suggested that these mutations may have effect white coat colors. However, the determination of white coat color followed the classical rule of epistatic effects with Agouti locus, especially wild type allele (A^wt) for this locus (or masked other unknown loci) [3,31].

As our previous hypothesis, the predominant black and red coat colors may carry E^D and e alleles at extension locus, respectively. However, the present study did not completely support the hypothesis, due to conserved putative alleles not being observed in all loci. All three missense mutations gave different frequency patterns for normal and mutated alleles, and this pattern resulted in no expression of animals carrying wildtype allele (E⁺) in investigated Mongolian goats (Supplementary Table S5). This was also likely due to animals carrying such allele were not included in this study. In some previous studies, the complete wild type phenotypes were described in Camosciate delle Alpi and Saanen for extension and agouti (A^wt) loci, respectively [3,32]. These wild type animals tended to carry a high frequency of normal alleles at the all loci. Fontanesi et al [3] investigated Derivata di Siria (red) goats in a previous study and the pheomelanin phenotype red, which they found and reported for clarification the dominance and incomplete penetrance of some “red alleles” at the extension locus.

Our findings as well as previous studies reported the highly conserved SNPs in wide range of goat populations may have an incomplete function on coat color phenotypes. This may be due to polygenic inheritance on coat color expression and masked effect of other loci. Therefore, previous studies [3,17, 31–34] emphasized that especially epistatic interaction of Extension and Agouti loci plays a significant role for this phenomenon. For an example, expression of Agouti alleles likely depends on at least one wild type copy (E⁺) putative allele at the Extenstion locus. More in detial, combination of dominant agouti (A^wd) and recessive extension (e) has been proposed for black phenotype [35].

In conclussion, we report SNPs in the MC1R gene of Mongolian goats in this study, and looked for their association with genotypes. Some of our findings on coat colors suggested that mutations in MC1R gene may have some role for determining coat color phenotypes. However, according to entire results obtained in this study also demonstrated that this locus is probably not the only factor for coat color trait. Further, the necessity of elucidate the genetic factors, which regulate the pigmentation of coat colors, especially interaction between the loci that associated with these traits were highlighted. Furthermore, the limitation for this study was a lack of sample from all registered breeds or strains of Mongolian goat, and this lack of samples may have resulted in the non-expression of previously reported and novel altered point mutations. Although, we believe that the largest number of the recorded coat colors of Mongolian goats were investigated in this study.