Functional divergence of the pigmentation gene melanocortin-1 receptor (MC1R) in six endemic Macaca species on Sulawesi Island

Coat color is often highly variable within and between animal taxa. Among hundreds of pigmentation-related genes, melanocortin-1 receptor (MC1R) plays key roles in regulating the synthesis of the dark eumelanin and the red–yellow pheomelanin. The six species of macaques that inhabit Sulawesi Island diverged rapidly from their common ancestor, M. nemestrina. Unlike most macaques, Sulawesi macaques commonly have a dark coat color, with divergence in shade and color pattern. To clarify the genetic and evolutionary basis for coat color in Sulawesi macaques, we investigated the MC1R sequences and functional properties, including basal cAMP production and α-MSH-induced activity in vitro. We found fixed non-synonymous substitutions in MC1R in each species. Furthermore, we found that six species-specific variants corresponded with variation in agonist-induced and basal activity of MC1R. Inconsistent with the dark coat color, four substitutions independently caused decreases in the basal activity of MC1R in M. hecki, M. nigra, M. tonkeana, and M. ochreata. Selective analysis suggested MC1R of M. nigra and M. nigrescens underwent purifying selection. Overall, our results suggest that fixed differences in MC1R resulted in different functional characteristics and might contribute to divergence in color among the six Sulawesi macaque species.


Results
Genotyping and selection analysis. We determined the MC1R nucleotide sequences to investigate polymorphisms in six endemic Macaca species on Sulawesi Island. From 51 individuals, we identified 15 alleles based on combination of 26 single-nucleotide polymorphisms (SNPs) (Supplementary Table S1), including 10 previously reported nonsynonymous SNPs 15 (Table 1). In particular, each species exhibited distinct fixed amino acids. All 10 nonsynonymous substitutions (P2R, P22L, M38V, G104S, H153P, M199L, C267Y, I293V, E304G, and R306C) were responsible for six species-specific haplotypes. Sulawesi macaques were separated into two clusters, one is the northern cluster, including M. nigra and M. nigrescens, the other one is the southern cluster, including M. hecki, M. tonkeana, M. maurus and M. ochreata (Fig. 2). Ancestral sequences of northern and southern cluster were constructed. Altogether, 9 of the 10 amino acid differences were species-specific, and one Each species showed specific nonsynonymous substitutions in a different part of MC1R. In M. hecki, all individuals shared two specific substitutions, C267Y and I293V, located in the third extracellular loop (EL3) and the seventh transmembrane region (TM7) of the receptor, respectively (Fig. 3). In M. nigra, the specific E304G substitution resulted in a change from a negatively charged residue to an uncharged residue in the C terminal domain of the receptor. In M. nigrescens, the P2R substitution resulted in a change from a nonpolar to positively charged residue in the N terminal region of the receptor. In M. tonkeana, the conservative substitution G104S resulted in a change from a nonpolar to polar residue in EL1. M. maurus possessed the most amino acid substitutions, including P22L, H153P, and M199L. SIFT (Sorting Intolerant From Tolerant; ≤ 0.05) and PROVEAN (Protein Variation Effect Analyzer; ≤ − 2.50) analyses consistently showed that five (G104S, H153P,  15 . To detect positive selection in MC1R in Sulawesi macaques, we applied three models using the PAML package: a branch model, site model, and branch-site model. First, using a site model to investigate selection throughout the phylogeny based on RADseq data 19 , we did not detect sites under positive selection, details seeing in Supplementary Table S3. We further calculated ω values for the lineages with melanism (M. nigra and M. nigrescens), the other four Sulawesi macaques, and M. nemestrina and M. mulatta using a branch model to examine selective constraints on the black coat color lineage. The ω values for the melanism lineage (ω = 0.086) were significantly lower than those for other Sulawesi macaques (ω = 0.968) and the lineage including M. nemestrina and M. mulatta (ω = 0.856) ( Supplementary Fig. S1). These results suggested that species with melanism underwent purifying selection. We further examined positively selected codon sites by a branch-site model. The null model was not rejected, and no positively selected sites were found.   www.nature.com/scientificreports/ Intracellular cAMP production under various concentrations of α-MSH is presented in Fig. 5. The responses of MC1R to agonist α-MSH differed substantially among species. α-MSH dose-dependently activated MC1R of all species, except for M. hecki MC1R. M. nigra MC1R showed a significantly lower maximal cAMP production than those for MC1R of other species (pairwise t-test, P < 0.05, BH-adjusted; Table 2). The  (Table 2). Exceptionally, M. hecki MC1R showed cAMP accumulation with 100 nM α-MSH stimulation. Because saturation was not reached, we could not determine the EC 50 values (Fig. 5c).
To determine the key residues affecting MC1R function, we designed several key mutants of MC1R based on SIFT and PROVEAN analyses. To understand the functional changes in both northern and southern clusters, we constructed vectors of the predicted ancestral northern (G304E mutant of M. nigra) and ancestral southern  Table 2. Summary of basal cAMP production and maximun cAMP production. 'Basal cAMP production' is the cAMP accumulation without agonist stimulation; 'Maximum cAMP production' is the cAMP accumulation which reaches saturation in response to 100 nM α-MSH. The results are presented as mean ± SEM obtained from at least 3 times independent experiments. ND indicates 'not determined' .

Discussion
We observed that the amino acid sequences of MC1R are conserved within each species (n = 10) and identified fixed species-specific amino acid substitutions in six closely related Macaca species. We further analyzed the functional features of MC1R itself and in response to the agonist α-MSH by a cAMP assay. We found that all six species-specific MC1R variants exhibited divergent basal activity and agonist-induced cAMP performance compared with those of the predicted ancestral sequences of the northern and southern clusters and M. nemestrina MC1R. And we identified the key residues responsible for MC1R function by site-directed mutagenesis.
We observed low nucleotide diversity in MC1R in each species of Sulawesi macaques, with an average of 0.067 × 10 -2 , which was similar to the estimate for human MC1R in African populations (π = 0.07 × 10 −2 ) 21 and three times lower than MC1R diversity in both M. nemestrina (π = 0.21 × 10 −2 ) and M. fascicularis (π = 0.20 × 10 −2 ; Supplemental Table S5). Our results suggested that the low MC1R diversity of M. nigra and M. nigrescens was a consequence of purifying selection for a dark coat color. Purifying selection on MC1R appears to be common in non-human primates 22 . Moreover, Nakayama and colleagues compared the nucleotide sequences of MC1R in 18 Macaca species and concluded that the gene was under purifying selection in the ancestral lineage of macaques and the silenus group 15 . In the present study, we verified that the amino acid sequence of MC1R is www.nature.com/scientificreports/ highly conserved in each Sulawesi macaque based on large population sizes and obtained additional functional information for each haplotype. Numerous non-synonymous substitutions in MC1R have been identified in mammals. In humans, many loss-of-function variants are associated with red hair color 23 . Notably, the R151C, R160W, and D294H mutations in MC1R are strongly associated with pale skin and red hair in Eurasians. Opposite to the dark coat color in Sulawesi macaques, we found that basal activity levels of MC1R were dramatically diminished in M. hecki, M. ochreata, M. nigra, and M. tonkeana, resulting from C267Y, R306C, E304G, and G104S mutations, respectively. Exceptionally, M. maurus has the most ancestral coat color (brownish) 16 . Consistent with this observation, M. maurus MC1R displayed comparable basal activity and agonist-induced activity to those of M. nemestrina MC1R. The increase of constitutive activity in M. maurus can be explained by a change in the secondary structure by the H153P substitution in IL2. Accordingly, the five species-specific mutations (G104S, H153P, C267Y, E304G, and R306C) in Sulawesi macaques identified in the present study independently affected the basal activity of MC1R.
The most influential substitution was C267Y (M. hecki) with respect to agonist binding activity. C267 and C275 form disulfide bonds between TM6 and TM7 affecting the tertiary structure of the receptor 24 . Single point mutation of C267 to glycine results in a complete loss of NDP-MSH binding. However, the serine mutant retained some agonist binding ability, though weaker than that of the wild-type locus 25 . Similarly, we found that C267Y exhibited agonist-independent basal activity and a cAMP response to α-MSH with an extremely high threshold (> 100 nM). In humans, the plasma concentration of α-MSH is 21.30 ± 0.63 nM 26 . Our results suggested that the natural C267 mutant with the hydrophobic residue tyrosine in M. hecki retained low agonist-induced activation; however, activity might be constant under the low concentration of α-MSH in the body.
MC1R EL1 is small compared with most GPCRs and deletions of EL1 resides are associated with melanism in the gray squirrel 27 , jaguar, and jaguarundi 28 . In the present study, the G104S variant (M. tonkeana) exhibited decreased basal activity and a slight right-shift in agonist-induced activity. This variant has also been detected in gibbons 22 and buffaloes 29 . Miao et al. speculated that the black coat color is associated with the allele carrying the G104S substitution 29 , but do they cause melanism is unclear. In silico functional prediction suggested that the G104S substitution is deleterious in buffaloes, and this was further supported by our functional results for the G104S variant in M. tonkeana. Our result suggests that G104S substitution could decrease MC1R activity and hence it is not likely to be causative for the melanism.
It is not clear how these species-specific substitutions in MC1R became fixed in Sulawesi macaques. MC1R acts as a genetic switch that determines whether dark eumelanin or light pheomelanin is produced for the regulation of coat color. Loss-of-function mutations usually lead to lighter, yellowish color, including in humans. In the present study, we found that most species-specific mutations resulted in reducing MC1R activity in vitro. Based on M. nemestrina as the extant ancestor, the predicted ancestral haplotypes of southern and northern clusters showed significantly reduced basal activity, with further reductions occurring independently in M. nigra, M. tonkeana, M. hecki, and M. ochreata (Supplementary Fig. S1). The basal activity and agonist-induced activity evolved independently in each Macaca species on Sulawesi Island, consistent with results for lemurs and platyrrhines 14 . Haitina et al. suggested catarrhine primates displayed conservation of dose-dependent α-MSH binding and activation, with variation in basal activity. However, we found that C267Y in M. hecki almost led to loss of α-MSH-induced cAMP production. The changes in MC1R function caused by these novel mutations are not simply limited to melanism in coat color in Sulawesi macaques, consistent with the results of Haitina et al. 14 and Nakayama et al. 15 . For example, we did not investigate the promoter region of MC1R or the sequences and expression levels of other pigmentation-related genes (e.g., POMC (proopiomelanocortin) and TYR (Tyrosinase)). Increased levels of α-MSH might lead to increased melanin production in the island population 30 . The limit of in vitro heterologous expression system using HEK293 cells might be one of the causes for the discrepancy between the protein molecular properties and pigmentation phenotypes in the monkeys. There might be a systematic difference in transfection efficiency, protein localization, dimerization, and internalization in vitro assay system with transient transfection. Further explorations of the expression levels and regulatory regions of MC1R and other pigmentation-related genes in vivo are needed to explain coat color divergence in the Sulawesi macaques.
Sulawesi macaques differ in patterns of unmelanized or light part on the forearm, cheek, upper arm, thighs and hind-shanks. In rhesus macaques, Bradley et al. did not find significant differences in the patterns of gene expression in comparing dark, intermediate and light hair 31 . Their results suggested the coat color variation from light to dark in rhesus macaques was unlikely to be due to differences in expression levels of key pigmentation genes; MITF (Melanocyte Inducing Transcription Factor), MC1R, MGRN1 (Mahogunin Ring Finger 1), ATRN (Attractin), SLC24A5 (Solute Carrier Family 24), TYRP1(Tyrosinase-related Protein 1) and DCT (Dopachrometautomerase) 31 . Hence, the expression of MC1R might be not causative to fine tuning of pattern difference in the case for Sulawesi macaques. ASIP gene plays a key developmental role in color patterning. Spatio-temporal regulation of ASIP can further modify MC1R activity. Mundy and Kelly suggested that mutations in ASIP coding region were not involved in color changes among closely related primate species 32 . Allele-specific expression of ASIP in body part has been found to be responsible for color pattern differences in mice 33 . We speculate that expression and regulatory differences at ASIP might play an important role in pattern variation in Sulawesi macaques. A protein expression analysis would further improve our understanding of variation among species and body parts.
Sulawesi macaques are morphologically differentiated, though individuals with intermediate traits, presumably hybrids, have been reported in the border zone of each species' distribution [34][35][36][37] . Gene flow between Sulawesi macaque species has presumed from the intermediate morphological characteristics 38 . The possibility of gene flow in a hybrid zone is also supported by an analysis of microsatellite markers 39 . Previous genetic analyses of hemoglobin 18 and TAS2R38 40 identified shared common haplotypes among Sulawesi macaques; however, we found that MC1R shows species-specific variants, without shared haplotypes. Divergence in MC1R was the Scientific Reports | (2022) 12:7593 | https://doi.org/10.1038/s41598-022-11681-z www.nature.com/scientificreports/ greatest between M. maurus and M. nigra, suggesting a high correlation with geographical distance. While the relationship between the protein sequence and coat color phenotype is still unclear, MC1R could be a speciesspecific marker gene for Sulawesi macaques. In addition to the apparent change-of-function MC1R variants for cAMP production in Sulawesi macaques, MC1R may have other important functions that are unique to each species.  Table S1) collected for a previous study 40 . The saliva was scrubbed with cotton swabs and samples were stored in 2 mL tubes with 1 mL of lysis buffer consisting of 0.5% sodium dodecyl sulfate, 100 mM ethylenediaminetetraacetic acid, 100 mM Tris-HCl, and 10 mM NaCl at room temperature 42  To amplify and sequence the entire coding region of the MC1R gene, primers (MC1R-F: 5′ ATG AGC TAA GCA GGA CAC C 3′; MC1R-R: 5′ CAA CAC CTT CAG AGG TCA GT 3′) were designed using the Primer3Plus website (website: http:// www. bioin forma tics. nl/ cgi-bin/ prime r3plus/ prime r3plus. cgi). MC1R was amplified using ExTaq DNA Polymerase (Takara Bio Inc., Shiga, Japan) by PCR under the following conditions: initial denaturation at 94 °C for 10 min, 45 cycles of denaturation at 94 °C for 10 s, annealing at 56 °C for 30 s, and extension at 72 °C for 1 min, followed by a final extension at 72 °C for 10 min. PCR products were sequenced using BigDye Terminator v. 3.1 (Applied Biosystems, Carlsbad, CA) and the sequencing products were separated by capillary electrophoresis using a 3130xl Genetic Analyzer (Applied Biosystems).

Reconstruction of MC1R haplotypes and selection analysis.
Sequences of intact MC1R of all the samples were aligned using MEGA X. A maximum-likelihood (ML) of the six species was reconstructed with 1000 bootstrap replicates. The ancestral amino acid sequences of each clusters were inferred using ML-based ancestral reconstruction in MEGA X. Multisite haplotypes were reconstructed from sequence data using DnaSP v. 6.12.03 43 . The genealogical relationships among haplotypes were constructed, rooted with the MC1R sequence of Papio hamadryas (Accession number AY205105.1), using the median-joining algorithm implemented in PopART 44 . To test for a signal of selection in the partial MC1R gene or lineage, we used the CODEML program in PAML 4 45 . The ratio ω (d N/d S) is a measure of selective pressure, where ω = 1, ω < 1 and ω > 1 correspond to neutral evolution, purifying and positive selection, respectively. Three models, branch, site and branch-site models were used to analyze the selection of MC1R in Sulawesi macaques. The phylogenetic tree was referred to RADseq data in Evans et al. 19 . Firstly, we used site models (M0: null, M1: nearly neutral selection, M2: positive selection, M3: discrete, M7: beta, and M8: beta and ω > 1) to determine candidates of positively selected sites. Model fit was evaluated using likelihood ratio tests (LRT). Subsequently, to test whether the ω ratio was different among lineages, the melanism lineage (M.  , and 10 μL of Opti-MEM per well. The transfection mixtures were preincubated for 20 min at room temperature before their addition to plates. Then, the plates were incubated at 37 °C under 5% CO 2 for 24 h. The cAMP assays were performed using the cAMP-Gs Dynamic Kit (Cisbio, Codolet, France). Briefly, on the day of the cAMP assay, α-MSH was diluted with Stimulation Buffer 1 to obtain a final concentration range of 10 −12 to 10 −7 M. Then, 7500 cells were added to each well of a white-walled 384-well plate and stimulated with different concentrations of α-MSH (Sigma-Aldrich, St. Louis, MO) for 30 min at room temperature. The cells were also stimulated with 20 µM forskolin (final concentration 10 µM; Cisbio) to normalize for cell numbers 47 . Following the protocol for the cAMP Dynamic Kit, fluorescence signals at 665 and 620 nm were detected using the FlexStation 3 Microplate Reader (Molecular Devices Japan, Inc., Tokyo, Japan). Values are expressed as ΔF (using ratios of 665 nm/620 nm for the assay and 665 nm/620 nm for the mock). Data are expressed as ΔF/ΔF max 48 , which is the ratio of the ligand-dependent increase (ΔF) to the maximal production of cAMP under activation by 20 µM forskolin (ΔF max ). ΔF/ΔF max values were fitted by the function f(x) = min + (max − min)/(1 + (x/EC 50 ) h )), where x is the ligand concentration and h is the Hill coefficient, using the drc package in R 49 . At least three independent measurements were conducted for each vector. Data are reported as mean values ± standard error of the mean (SEM).

Data availability
DNA sequences are available in DDBJ; DDBJ accessions LC632229 to LC632296.