Functional colour genes and signals of selection in colour‐polymorphic salamanders

Coloration has been associated with multiple biologically relevant traits that drive adaptation and diversification in many taxa. However, despite the great diversity of colour patterns present in amphibians the underlying molecular basis is largely unknown. Here, we use insight from a highly colour‐variable lineage of the European fire salamander (Salamandra salamandra bernardezi) to identify functional associations with striking variation in colour morph and pattern. The three focal colour morphs—ancestral black‐yellow striped, fully yellow and fully brown—differed in pattern, visible coloration and cellular composition. From population genomic analyses of up to 4,702 loci, we found no correlations of neutral population genetic structure with colour morph. However, we identified 21 loci with genotype–phenotype associations, several of which relate to known colour genes. Furthermore, we inferred response to selection at up to 142 loci between the colour morphs, again including several that relate to coloration genes. By transcriptomic analysis across all different combinations, we found 196 differentially expressed genes between yellow, brown and black skin, 63 of which are candidate genes involved in animal coloration. The concordance across different statistical approaches and ‘omic data sets provide several lines of evidence for loci linked to functional differences between colour morphs, including TYR, CAMK1 and PMEL. We found little association between colour morph and the metabolomic profile of its toxic compounds from the skin secretions. Our research suggests that current ecological and evolutionary hypotheses for the origins and maintenance of these striking colour morphs may need to be revisited.


| INTRODUC TI ON
Animal coloration is conspicuously affected by natural and sexual selection (Caro, 2005) and therefore is an ideal system in which to study the evolutionary processes that generate and maintain biological diversity (Andrade et al., 2019;Barrett et al., 2019;Harris et al., 2020;Hubbard, Uy, Hauber, Hoekstra, & Safran, 2010;Kusche, Elmer, & Meyer, 2015;Stevens & Ruxton, 2012). Within the terrestrial vertebrates, amphibians present some striking colour patterns that vary both within and across species and are often presumed to be associated with ecological and physiological functions (Cott, 1940;Thayer, 1909). For example, putative antipredator strategies that are based on colour include crypsis or disruptive coloration (concealment through background colour matching or camouflage by pattern) and aposematism (an interspecies warning coloration communicating toxicity or unpalatability; Kerr, 1919;Ruxton, Sherratt, & Speed, 2004). Colour variation may also be associated with thermoregulation and water loss reduction (Duellman & Trueb, 1994;Moore & Ouellet, 2015;Rudh & Qvarnström, 2013), for example as suggested by an observed increase in melanic morphs at high altitudes and latitudes (Alho et al., 2010;Vences et al., 2002Vences et al., , 2014. However, inferences about the selective pressures that generate and maintain colour polymorphisms in amphibians have often been speculative and correlational, and the molecular underpinnings poorly understood (Hoffman & Blouin, 2000;Rudh & Qvarnström, 2013).
The ecological and evolutionary genomics of amphibians has been somewhat impeded by logistical and biological complications compared to the rapid recent advances in mammals, fish and birds (McCartney-Melstad, Mount, & Shaffer, 2016;Rudh & Qvarnström, 2013). For example, the challenges of captive breeding hamper genetic mapping for quantitative traits, especially in those species that have internal fertilization, long-term sperm storage, lengthy gestation times and small clutch sizes (Rudh & Qvarnström, 2013). Second, in many cases the colour variations are among species rather than within, which hinders distinguishing genomic patterns associated with colour rather than species differences. Third, amphibians have very large genome sizes, ranging up to 121 pg (Gregory, 2019), which has slowed the pace of reference genome development and functional genomic analyses compared to other vertebrates (Koepfli, Paten, & O'Brien, 2015). Large genome size in salamanders is due to an excess of repetitive elements but the size and number of coding genes is similar to that of other vertebrates (Nowoshilow et al., 2018). To advance knowledge of how selection affects the origins and maintenance of amphibian colour variation, analyses of colour polymorphisms and related genetic patterns in environmental context are needed.
Here, we provide novel insights into the genetic basis of coloration in amphibians using a highly colour-polymorphic lineage of salamander. Fire salamanders (subspecies in Salamandra salamandra) have an array of yellow, black, brown, red, orange, white, striped and spotted colorations that are usually fixed within subspecies (Sparreboom & Arntzen, 2014;Thiesmeier, 2004;Thorn & Raffaëlli, 2001). The striking colour patterns are widely assumed to be aposematic, functioning as warning coloration (Beukema, Speybroeck, et al., 2016). Fire salamanders are able to release endogenously produced noxious secretions from their enlarged paratoid glands at the base of the head and smaller glands along the spine as an antipredator defence, which are quite toxic compared to those of most European amphibians (Lüddecke, Schulz, Steinfartz, & Vences, 2018;Trochet et al., 2014;Vences et al., 2014). Alternatively but nonexclusively, the black and bright spotted or striped coloration may act as a form of concealment in forested habitats, as aposematic coloration can also be visually disruptive (Stuart-Fox & Moussalli, 2009), although this has yet to be empirically tested. Some Salamandra lineages have evolved melanic colorations by losing the yellow components (Bonato & Steinfartz, 2005;Steinfartz, Veith, & Tautz, 2000;Vences et al., 2014). This is presumed to be due to selection for thermoregulation because, being ectotherms, darker skin at higher altitudes protects against ultraviolet light and enables more efficient warming from the sun. Repeated independent evolution of melanic and striped morphs and sister lineages at higher altitude are suggestive evidence for this hypothesis (Beukema, Speybroeck, et al., 2016;Lüddecke et al., 2018;Rodríguez et al., 2017;Vences et al., 2014). Maintenance of colour morphs within a population might be promoted by colour-based assortative mating, as proposed in other salamanders (Acord, Anthony, & Hickerson, 2013), but not yet tested in this species. Aposematism, crypsis, thermoregulation and assortative mating have been demonstrated to different extents in colour polymorphism generation and maintenance at the phenotypic level but to date have been rarely evaluated in a salamander.
In this study we aimed to identify if signals of response to natural selection are associated with coloration in colour-and pattern-variable salamanders. Experimentally we used a lineage of fire salamanders that have striking solid-coloured yellow and brown colour morphs and intermediate colorations alongside the more typical ancestral yellow-black striped morph (Beukema, Nicieza, Lourenço, & Velo-Antón, 2016;Köhler & Steinfartz, 2006; Figure 1). Our expectation is that functional genomic variation differs between colour-polymorphic phenotypes and there are genomic signals of response to selection. We: (a) characterized coloration for pattern, visible and nonvisible colour, and cellular K E Y W O R D S amphibian, coloration, colour candidate genes, gene expression, genetic architecture, genomic association, population genomics, response to selection, toxicity structures; (b) associated genome-wide genetic variation with colour morph and inferred signals of selection on genomic patterns; (c) identified functional candidate genes for coloration based on expression; and (d) assessed the intersections between genome-wide association, functional expression and signals of selection. Using this, we evaluated support for three major standing hypotheses for the existence of colour morphs: local adaptation to altitudinal differences, assortative mating and differential antipredator toxins (Acord et al., 2013;Beukema, Speybroeck, et al., 2016). The coexistence of different colorations within a single environment offers a natural experiment for inferring genetic bases, ecological and phenotypic correlates of colour, and disentangling response to selection.

| Field sampling
Adult salamanders were sampled at seven localities in Asturias, Spain, from 78 to 1,312 m asl ( Figure 1). Striped individuals were found at all sites; brown and yellow salamanders were found at localities 4-7. Gradations between striped, yellow and brown also exist in the population but our sampling focused on characterizing the most representative colorations. Standardized dorsal images were taken using a Canon EOS 70D DSLR camera using a 100-mm lens and with an X-Rite ColorChecker. Tissue samples (toe or tail clips) were stored in ethanol. Salamanders were released where captured, except those representatives sampled for microscopy and gene expression studies. Sex was inferred from cloacal swelling in males.

| Colour pattern analysis
While salamanders can be clearly grouped into colour morph by eye, colour patterns were additionally examined objectively and quantitatively from a single densely sampled colour polymorphic sample site (Site 7) to mitigate any interpopulation effects. Photographs were standardized in Adobe Photoshop CS6: (a) a custom camera profile (generated using X-Rite's DNG ProfileManager) was applied, (b) the white balance of each image was set using the white square of the X-Rite ColorChecker, and (c) the exposure was modified until the RGB (red, green and blue) values of the white ColorChecker square equalled 240:240:240 (±5).
Dorsal colour patterns were analysed in r (R Core Team, 2013) using Patternize (Van Belleghem et al., 2018). Only clear images without obvious body contortions were used: 61 striped, 115 yellow and 44 brown (total N = 220). The ratio of females to males per colour morph was 1.3, 1.3 and 1.6 respectively. Image registration and colour pattern extraction was carried out using the patRegK function. The patPCA function was used to conduct principal component analyses (PCAs) to look at differences between the three colour morphs. Resulting PC scores were compared via one-way ANOVA.

| Transmission electron microscopy (TEM)
To identify the chromatophores that underlie different skin colours, two dorsal skin samples were dissected from striped, yellow and brown salamanders ( Figure S1b; N = 2 individuals per colour morph).

| Metabolomics of toxin
We investigated the relationship between colour phenotype and the metabolomic content of the toxic secretions using gas chromatogra-
Loci were assembled de novo in stacks version 1.44, with the minimum reads required to create a stack as six (other settings left on default). Genotype and haplotype corrections in individual samples were then conducted using the stacks rxstacks pipeline: loci for which 25% or more individuals had a confounded match to the catalogue were removed, excess haplotypes were pruned, sequencing errors were removed using the bounded model and single nucleotide polymorphisms (SNPs) were recalled, and catalogue loci with an average log likelihood less than -10.0 were removed (Rochette & Catchen, 2017). Following this, the catalogue of loci was rebuilt and each sample was mapped against it. The export_sql.pl function was used to generate a whitelist of loci containing one or two SNPs. Samples were filtered through the stacks populations pipeline as a single population, with whitelisted loci retained if they were present in ≥75% of samples, had a minimum in-
The best-fit value was determined using structure harvester (Earl & VonHoldt, 2012) to identify the highest ΔK (following Evanno, Regnaut, & Goudet, 2005). F ST among sample sites and between colour morphs in sample sites 6 and 7 was calculated through pairwise differentiation tests in genodive (Meirmans & van Tienderen, 2004); a post-hoc Bonferroni correction for multiple testing was applied to p-values using the p.adjust function in the r stats package (R Core Team, 2013).

| Gene expression
We used transcriptomic analyses to identify gene expression in yel- Raw read pairs ranged from 20.17 to 29.94 million per sample, with 86%-95% retained after quality filtering in trimmomatic version 0.36 (Bolger, Lohse, & Usadel, 2014

| Colour genomics
To reduce the impact of population structure on the inference of colour-associated loci across sampling sites, data for the 80 genotyped individuals used for population genetic analyses were refiltered through stacks populations as above, but without missing data constraints. We excluded sites deviating from Hardy-Weinberg equilibrium and containing more than 25% missing data across all individuals using vcftools, generating a data set of 2,149 SNP loci.
Missing data (13.4%) were imputed using fastphase version 1.4.8 (Scheet & Stephens, 2006) and SNP loci were converted into binary format using plink (Purcell et al., 2007). Genotypes were then corrected for population structure through linear regression with the first 10 PC scores from a PCA of the data in the R package adegenet.
Random forest (RF) analyses were conducted in r using the randomForest package (Liaw & Wiener, 2002). Briefly, three independent analyses of 100,000 trees were run (ensuring a Pearson correlation coefficient of ≥0.95 between runs), which provided a mean permuted importance statistic (%IncMSE) for each SNP locus. The top 101 loci-roughly corresponding to the "upper end of the elbow" of the importance distribution, following Laporte et al. (2016)-were then subset and a backwards purging process used to determine the cohort that explained the highest amount of variance in the data.
Briefly, three analyses of 10,000 trees were run, the mean r 2 across runs was calculated and the SNP with the smallest %IncMSE was removed. This was repeated until the data set contained only two SNP loci, with the iteration displaying the highest r 2 indicating the set of covarying loci best able to discriminate between colour morphs. Latent factor mixed model (LFMM; Frichot, Schoville, Bouchard, & François, 2013) analyses were conducted in the r package LEA (Frichot & François, 2015). This package tests associations between loci and environmental (or phenotypic) gradients using an MCMC algorithm for regression analysis. The analysis was run 100 times, increasing the power to detect true associations, assuming one latent factor (K; an exact value is not essential as it is estimated). Z-scores (the number of standard deviations above or below the population mean a data point is) from these runs were combined and the median value per locus calculated. The genomic inflation factor (λ) was then

| Colour morph differences
Objective quantification of colour pattern supported the grouping of individuals into the three colour categories: black-yellow striped, yellow morph, or brown morph (Patternize PC1 scores: F (2, 217) = 429.8; p < .005; Figure 2a). Colour patterns did not differ between the sexes (F (1, 218) = 0.011; p = .916) so data from males and females were combined. The three colour morphs differed in hue at the head, paratoid, dorsal and lateral landmarks (Table S1). Colour reflectance varied clearly across skin that was brown, yellow or black ( Figure 2b); for example, the yellow on black-yellow striped individuals is the same yellow as on the yellow morph, and brown on the brown morph is similar to black on striped individuals (Table S1; Figure S2).
The types and densities of chromatophore cells in skin sections from dermal and epidermal layers at different colour landmarks were qualitatively different ( Figure 2c). Based on visual inference, yellow skin had xanthophores, which are yellow-pigment-containing, in the epidermal layer and a high density of light-reflecting iridophores in the dermal layer. In poikilothermic skin, xanthophores contain yellow pigment while iridophores are structurally light-reflecting (Bagnara, 1972;Bagnara & Ferris, 1971;Bagnara, Frost, & Matsumoto, 1978;Kimura et al., 2014). In contrast, in both skin layers of brown and black skin only melanophores were evident, which are dark pigment-containing cells (Bagnara, 1972;Bagnara & Ferris, 1971;Bagnara et al., 1978;Mills & Patterson, 2009). The density of melanophores in the epidermal and the dermal layers was higher in the black than brown skin. In yellow skin, melanophores were infrequent and found as single cells or localized dermal clusters. These results demonstrate that yellow and black/brown skin colours are due to different underlying cellular structures.

| Population genomics
We tested for genome-wide population genetic differentiation between the different colour morphs based on 4,702 loci. We found low but significant genetic differentiation between geographical localities regardless of colour (F ST : 0.018-0.196, Table S3). Within the colour-polymorphic sample localities, there was no neutral genetic differentiation between colour morphs (all colour morph comparisons: F ST < 0.003, p > .12; Table S2). Population clustering reflecting genomic proportions was driven by sample localities and not colour ( Figure 3b). Individuals of different colours are intermixed in a maximum-likelihood (ML) tree and group into subclades by population ( Figure 3a), as also evidenced in principal components ( Figure S3).
Overall these suggest no structuring by colour morph.

| Genotype-phenotype associations and functions
We tested for an association between genomic variation and phenotype and identified 21 loci best able to discriminate between individuals based on their colour morph using an RF analysis (61.22% variance explained; 2,149 loci in the data set; Figure S4a; Table 1).
Five of these could be mapped indirectly via salamander contigs to vertebrate genes (Table S4) Table S5). Of these, 5-9 loci overlapped between pairwise comparisons (Table S5). No F I G U R E 2 Colour phenotype characterization. (a) Patternize PC1 scores for salamanders after a priori assignment to colour morph (n = 220). (b) Average reflectance spectra of yellow (n = 14), brown (n = 9) and black (n = 17) skin across individuals (shading = standard deviation). (c) TEM images of Salamandra skin (circles in the bottom-right corner represent the skin colour): (1) TEM image of epidermal yellow skin from a striped salamander;

| Signals of selection
(2) epidermal brown skin from a brown salamander; (3) epidermal black skin from a striped salamander; (4) dermal yellow skin from a striped salamander; (5) dermal brown skin from a brown salamander;  (Table S5).
Of all loci showing a signal of response to selection, 10 overlapped with those identified as being associated with coloration through genotype-phenotype association analyses (Table S5).

| Skin colour gene expression
Based on gene expression differences between skin of different colour, we identified a total of 167 DE genes (adjusted  (Table S6; Figure S7). The striped individuals were sampled for yellow and black skin (one landmark each; Figure S1) and are within the latter comparison; those minimize confounding factors associated with interindividual variation and allow definitive identification of genes differentially expressed between colours.
Of all transcripts that were differentially expressed, 60 could be putatively associated with genes with known involvement in animal pigmentation through assessment by imp, amigo 2 and literature searching (Table S7): 44 with melanophores, five with xanthophores, four with iridophores and 12 with other pigmentation-related processes (with some overlap, see Table S7 for full details).

(a)
Based on over-representation analysis for genes up-regulated in skin of different colours, we identified molecular pathways putatively contributing to differences in expression (Table S8).
Genes up-regulated in both brown and black skin compared to yellow skin included GO terms related to pigmentation, melanin biosynthesis and melanosomes. Broadly, genes up-regulated in yellow skin compared to black were associated with inosine and urate metabolic processes. None of the GO terms identified between brown and black skin have been associated previously with animal pigmentation. This suggests that the expression of brown and black skin may result from as yet uncharacterized molecular and cellular factors.

| Toxic secretions and colour polymorphism
We identified 18 metabolites in toxic secretions shared across most or all individuals (Table S9). A PCA of metabolite intensities for S. s. bernardezi, which indicates the relative abundance and variation, showed that the profiles of striped, yellow and brown salamanders overlapped along the first axis (F (2,15) = 0.219; p = .806) but significantly separated along the second axis (F (2,15) = 5.595; p = .015; Figure 6a).
No significant differences were found in intensity between the colour morphs at either those four characterized alkaloids or at any of the 14 uncharacterized metabolites (Table S10; Figure 6b). There was a difference between the S. s. bernardezi (i.e., all three colour morphs combined) and the outgroup comparison S. s. terrestris at cycloneosamandione (Peak02_375), samandarone (Peak07_375) and samandarine (Peak06_449) (Table S10), which demonstrates the method is sensitive to variation and that these subspecies differed.
Among the three colour morphs, some trends were evident, such as lower levels of samandarine in the yellow morph (one-way ANOVA: F (2,15) = 3.11; p = .074). However, a large amount of variability was seen within and between the colour morphs ( Figure 6b).

| Molecular basis of colour in salamanders
Here we identified genomic loci associated with colour morph in these polymorphic salamanders that are yellow, brown or black-yellow striped. A total of 29 loci were inferred to covary with colour morph (21 from RF approaches and nine from LFMMs)-one locus overlapped across these complementary analyses (locus ID 28130 indirectly annotated to CAMK1). It is particularly striking that we inferred shared candidates across lines of analysis with transcriptomic and genomic data sets, and given the very large 34-Gb genome size of this species (Gregory, 2019), this suggests substantial regions in genetic linkage may underpin the genetic basis of colour in this species.
Although this is a nonmodel species with no reference genome, eight of those 29 colour-associated loci could be indirectly mapped to characterized genes by comparison with the biomedical model axolotl and other vertebrates. This included the colour gene TYR, which in mammals encodes the enzyme tyrosinase involved in the conversion of tyrosine into melanin pigments (Murisier & Beermann, 2006). Mutations in TYR have been shown to result in an albino phenotype in laboratory-reared axolotl salamanders (Woodcock et al., 2017), where melanophores develop but remain unmelanized (Frost, Epp, & Robinson, 1986). TYR is a particularly well-supported TA B L E 1 List of de novo loci that random forest and LFMM analyses inferred to be associated with colour morph, and putative gene IDs identified using the sal-site assembly version 4.0 with NCBI blast candidate from several lines of evidence in our study; yellow skin also showed significantly reduced TYR expression in the transcriptome analyses and very few melanophores were seen in yellow skin through microscopy.

Analysis
Several other candidates found to be significantly associated with colour morph also have biologically relevant links to coloration. For example, BAZ2A is known to share enhancers (GH12G055965 and GH12G056040) with the characterized colour gene PEML (premelanosome protein; Rebhan, Chalifa-Caspi, Prilusky, & Lancet, 1997).
The gene CAMK1 is hypothesized to be involved in guanosine-related processes in humans, with guanosine being a key constituent of iridophore cells (Ide & Hama, 1972). We found iridophore cells abundantly in the dermal layer of yellow Salamandra skin (Figure 2c). Revealingly, CAMK1 was found to be associated with colour morph in the RF and the linear model analysis, suggesting it warrants further investigation for its role in this colour diversity. However, overall these candidates are somewhat tentative because of the lack of a closely related reference genome. Future work is needed to validate the loci for coloration.
Our findings from gene expression are more direct associations of colour morph and molecular bases and these supported several aspects of the genomics inference. Notably, the four most down-regulated transcripts in yellow skin were related to three genes involved in melanin production (PMEL, TRYP1 and TYR). Melaninrelated genes were also important in distinguishing between yellow and brown skin, with the most up-regulated (DEBF family and TFEC) and down-regulated (MAPK12, PMEL, MLANA, CCM2L and PDE1B) genes in yellow skin all being melanin-process-related. Across vertebrates melanin is well established to be the main pigment of darker skin, fur and feather (Mills & Patterson, 2009). This strongly implicates melanin in the differences between yellow and dark (brown or black) skin in salamanders, consistent with our structural observations that melanophores are less abundant in yellow dermis (Figure 2c-4). The influence of melanin-pathway down-regulation compared to structural melanophore abundance in yellow skin warrants further study.
The seven transcripts most up-regulated in yellow skin compared to black skin mapped to genes with roles in pigmentation in vertebrates, including iridophore (PNP and ADA2), xanthophore (PAX7 and SLC2A11) and leucophore (SLC2A9) production. Cells of the dermal chromatophore unit found in amphibian skin (Bagnara & Hadley, 1973) are either pigment-containing (yellow-red xanthophores) or (2) 0 9  (Table S7). While the exact mechanisms remain unclear, it is known that keratinocytes interact with melanosomes either through melanin pigment transfer or uptake (Aspengren, Hedberg, & Wallin, 2006). Brown and black melanin pigments also share a common molecular basis, both being eumelanin polymers (Bagnara et al., 1978;Ito & Wakamatsu, 2011;Meredith & Sarna, 2006). Given this, we speculate that brown skin arises as a result of modifications in one or more melanin pathways, which either change the kinds or ratios of eumelanin pigments being produced in melanosomes. This would be consistent with recent studies on human skin pigmentation (Crawford et al., 2017;Ebanks et al., 2011) and could explain the cellular similarity, yet visual difference, between black and brown skin. Alternatively, the size and arrangement of melanophores in brown skin may be different to that in black skin, which could be elucidated by further quantitative studies on cellular structure.
In addition to specific genes, we found that higher level molecular pathways differed between the colour morphs. When compared to yellow skin, genes up-regulated in both brown and black skin included pigmentation-associated GO terms relating to melanin biosynthesis and melanosomes. Genes up-regulated in yellow skin compared to black were broadly associated with inosine and urate metabolic processes. This is functionally relevant because inosine is important in the production of xanthophores and iridophores (Krauss et al., 2013) and uric acid is important in the production of leucophores (Kimura et al., 2014). Our results suggest valuable insights into the biochemical pathways involved in chromatophore production, which remain poorly characterized (McLean, Lutz, Rankin, Stuart-Fox, & Moussalli, 2017) despite sharing some common developmental pathways with melanophores (see Beirl, Linbo, Cobb, & Cooper, 2014;Darias et al., 2013;Woodcock et al., 2017).

| Evidence colour is under selection?
Our findings from selection analyses give tentative insights into the mechanisms generating and maintaining this variation in colour and pattern. We identified 142 loci showing putative signals of response to selection between colour morphs. Several were shared across multiple pair-wise comparisons. Of the loci showing signals of response to selection, 10 overlapped with those associated with coloration through genotype-phenotype analyses (see Table S5).
Such overlap lends support to the inference (Axelsson et al., 2013) and suggests a functional link between colour loci and targets of selection.
Four loci that are candidates for response to selection are particularly noteworthy because they are functionally related to colour in other organisms. First, the CAMK1 locus showed decreased genetic diversity in both black-yellow striped and yellow salamanders compared to brown salamanders, suggesting that divergent selection is driving the differentiation between colour morphs in this locus. This is striking because based on our gene expression analysis we hypothesized that CAMK1 is involved in iridophore production. CAMK1 is involved in guanosine-related processes and guanosine is a key constituent in iridophores in frogs (Ide & Hama, 1972). These cells are found in yellow skin and apparently have been lost in brown individuals ( Figure 2). As reduced genetic diversity in one morph compared to another can be indicative of a selective sweep (Wilson Sayres, Lohmueller, & Nielsen, 2014), this suggests that CAMK1 may be under selection in some colour morphs. Furthermore, CAMK1 was found in multiple association and selection analyses in this study, demonstrating a robust and consistent signal. These data suggest a functional link between the signals of selection we inferred and structural pigmentation, although the target and mechanism need to be resolved.
The TYR locus was also found to be putatively under selection between striped and yellow salamanders. It showed reduced genetic diversity in the derived yellow and brown morphs, suggesting that it is involved in colour pattern changes.  (Mallarino et al., 2016) and grey versus black feather colorations in crows (Poelstra et al., 2014). Both of the derived salamander colour morphs-yellow and brown-are solid coloured and our findings suggest a possible role of TYR in the loss of stripes that warrants further testing.

| Ecological and evolutionary mechanisms of colour polymorphism
Post-hoc we can use the findings on phenotypes and genotypes to explore some standing hypotheses for the ecological and evolutionary drivers of amphibian colour pattern evolution: local adaptation to altitudinal differences, assortative mating or differential antipredator toxins (Beukema, Nicieza, et al., 2016;Beukema, Speybroeck, et al., 2016;Cott, 1940;Vences et al., 2014).
First, we do not find support for the hypothesis that the derived skin colour morph is adaptive for thermoregulation requirements at different altitudes, as had been suggested in other salamander lineages (Vences et al., 2014). We sampled salamanders from near sea level up to 1,312 m in the mountains and both the brown and the yellow colour morphs were always found alongside each other and were found at low to intermediate altitudes (172-679 m asl).
In several locations all colour morphs that we studied were found together ( Figure 1). This is in agreement with research spanning a broader distribution of S. s bernardezi and niche breadths (Beukema, Nicieza, et al., 2016). Our finding suggests that the derived colour morphs are not differentially found at high or low altitudes at the spatial scale of this study, although historical effects and processes for altitudinal adaptation operating at the origin of different colour morphs cannot be ruled out.
Second, we do not find support for the hypothesis that there is contemporary reproductive isolation between colour morphs.
Using comprehensive genome-wide markers we found that, within any colour polymorphic sample locality, there was no genetic differentiation between colour morphs. The lack of population genetic differentiation between colour morphs at genomic loci is in contrast to predictions; it had been suggested that these rare colour variants had a distinct history and previously were considered a separate subspecies (Beukema, Nicieza, et al., 2016;García-París, Alcobendas, Buckley, & Wake, 2003;Steinfartz et al., 2000). Our findings are also consistent with these derived yellow and brown morphs having arisen in situ within the ancestral black-yellow striped variant but without colour assortative mating. However, allopatric divergence followed by secondary contact with high gene flow between morphs cannot be ruled out.
Future studies with colour-informed demographic modelling and outgroup sampling are required.
Third, we find equivocal support for a relationship between colour morph and the chemical profile of the famous toxic skin secretions of fire salamanders (Lüddecke et al., 2018;Mebs & Pogoda, 2005). Some differences are suggestive, such as in samandarine, and warrant further testing. Studies on other species and subspecies of Salamandra found no correlation between toxic secretions and colour pattern (Preißler et al., 2019;Sanchez et al., 2019;Vences et al., 2014) but these studies did not include the yellow or brown colour morphs. Colour in salamanders may therefore differ from amphibian systems where levels of conspicuousness and toxicity have been shown to be tightly correlated (Amézquita et al., 2017;Maan & Cummings, 2012). However, the sample size used for metabolomic analyses in the present study is limited and further research on toxicity and aposematism is needed. Our investigations suggest we may need to consider new hypotheses for the generation and maintenance of colour variability in fire salamanders. The signals of selection found between these colour morphs do indicate a potentially adaptive role but in response to what selective pressure is not known. Definitively inferring the location and number of structural and regulatory changes underpinning the differences between salamander skin colour generation will require whole genome data and a well-annotated reference genome, something that is not currently available for Salamandra but will be valuable for more comprehensive testing.

| CON CLUS ION
Our study identifies several molecular associations with coloration in polymorphic salamanders and suggests that these rare colour morphs are under selection. Natural admixture across colour morphs enabled us to analyse the association of population structure and coloration, without confounding background structure.
The inferred genotype-phenotype associations, combined with an overlap of loci across different statistical approaches and genomic and transcriptomic data sets, enhances our confidence that we have in fact identified loci related to functional differences between colour morphs. These several lines of evidence from different 'omic data sets to known pigmentation genes, including TYR, CAMK1 and PMEL, implicate a role in different colorations and suggest important candidate genes for future studies. We suggest that the current hypotheses are not sufficient to explain the ecological drivers of this distinctive and geographically restricted colour diversity. This research will facilitate studies on the evolution and convergence of colour patterns across vertebrates, which in turn will allow us to assess the consistency-as well as constraints onvertebrate colour pattern evolution by natural selection.