
==== Front
Nat Commun
Nat Commun
Nature Communications
2041-1723
Nature Publishing Group UK London

39227572
51927
10.1038/s41467-024-51927-0
Article
PMEL is involved in snake colour pattern transition from blotches to stripes
http://orcid.org/0000-0002-0831-5221
Tzika Athanasia C. athanasia.tzika@unige.ch

1
http://orcid.org/0000-0002-8595-7703
Ullate-Agote Asier 13
Helleboid Pierre-Yves 1
Kummrow Maya 2
1 https://ror.org/01swzsf04 grid.8591.5 0000 0001 2175 2154 Laboratory of Artificial & Natural Evolution (LANE), Department of Genetics & Evolution, University of Geneva, Geneva, Switzerland
2 https://ror.org/02crff812 grid.7400.3 0000 0004 1937 0650 Tierspital, University of Zurich, Zurich, Switzerland
3 grid.508840.1 0000 0004 7662 6114 Present Address: Biomedical Engineering Program, Center for Applied Medical Research (CIMA), Universidad de Navarra, Instituto de Investigación Sanitaria de Navarra (IdiSNA), Pamplona, Spain
3 9 2024
3 9 2024
2024
15 765515 4 2024
20 8 2024
© The Author(s) 2024
2024
https://creativecommons.org/licenses/by/4.0/ Open Access This article is licensed under a Creative Commons Attribution 4.0 International License, which permits use, sharing, adaptation, distribution and reproduction in any medium or format, as long as you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons licence, and indicate if changes were made. The images or other third party material in this article are included in the article’s Creative Commons licence, unless indicated otherwise in a credit line to the material. If material is not included in the article’s Creative Commons licence and your intended use is not permitted by statutory regulation or exceeds the permitted use, you will need to obtain permission directly from the copyright holder. To view a copy of this licence, visit http://creativecommons.org/licenses/by/4.0/.
Corn snakes are emerging models for animal colouration studies. Here, we focus on the Terrazzo morph, whose skin pattern is characterized by stripes rather than blotches. Using genome mapping, we discover a disruptive mutation in the coding region of the Premelanosome protein (PMEL) gene. Our transcriptomic analyses reveal that PMEL expression is significantly downregulated in Terrazzo embryonic tissues. We produce corn snake PMEL knockouts, which present a comparable colouration phenotype to Terrazzo and the subcellular structure of their melanosomes and xanthosomes is also similarly impacted. Our single-cell expression analyses of wild-type embryonic dorsal skin demonstrate that all chromatophore progenitors express PMEL at varying levels. Finally, we show that in wild-type embryos PMEL-expressing cells are initially uniformly spread before forming aggregates and eventually blotches, as seen in the adults. In Terrazzo embryos, the aggregates fail to form. Our results provide insights into the mechanisms governing colouration patterning in reptiles.

Corn snakes are a useful model organism for studying skin pigmentation. Here they characterize the Terrazzo morph and identify the role of PMEL in the patterning of pigmented cells in the skin, explaining how the wild-type blotches change to stripes in Terrazzo.

Subject terms

Ectoderm
Evolutionary developmental biology
Herpetology
Development
https://doi.org/10.13039/501100000854 Human Frontier Science Program (HFSP) RGP0037/2022 Tzika Athanasia C. Swiss National Sciences Foundation - 310030_204466 Ernst and Lucie Schmidheiny Foundation - 10_2023 Fonds Général de l’Université de Genève - 23_28 Emile Plantamour Fund - 2024/12 iGE3 PhD award - 2015 Sara Borell Grant - CD22/00027issue-copyright-statement© Springer Nature Limited 2024
==== Body
pmcIntroduction

Blotches, spots, labyrinths, stripes, and any combination of these are among the patterns encountered in animal skin colouration. In-depth work has been done on classical model species to identify the key players involved in the establishment of these spectacular phenotypes1 and to understand how these elements interact2,3. In zebrafish, a species whose skin colouration has been extensively studied, it is the interactions of three types of cells in the skin—the melanophores, the xanthophores, and the iridophores—that produce the yellow and black stripes of the adult4. These cells originate from the neural crest and form the larva skin pattern during embryogenesis, and then the adult pattern during metamorphosis. The recent development of single-cell transcriptomics makes it possible to investigate the gene expression profile of each chromatophore type at different developmental stages, thus furthering our understanding of how they differentiate and mature5,6. By elucidating the mechanisms that bring about the wild-type zebrafish colouration, it is also possible to understand the phenotypic diversity of other teleosts7,8, but it remains to be seen if the extensive findings in this vertebrate lineage can be extrapolated to the other vertebrate lineages9.

Reptiles, comprising more than 12,000 species10, exhibit equally spectacular colouration phenotypes as teleost fish but the inherent difficulties in maintaining and breeding them in captivity, mainly due to their larger size and the long generation time, has restricted the number of studies on these animals. Yet a handful of squamate species (snakes and lizards) are currently promoted as new models. In recent years, several colour morphs have been characterised for the corn snakes (Pantherophis guttatus)11,12, the leopard geckos (Eublepharis macularius)13,14, and the ball pythons (Python regius)15,16. These animals have been kept in captivity for decades by private breeders and numerous spontaneously-occurring mutations have appeared, affecting both their colouration and their colouration pattern. The systematic genetic characterisation of these colour morphs, coupled with developmental, transcriptomic, and functional analyses, will help us discover the mechanisms that establish the patterning of squamate chromatophores.

In this work, we focus on Terrazzo, a corn snake colour morph that presents stripes instead of blotches (Fig. 1). We genetically characterize this morph using a mapping-by-sequencing approach. We shortlist 35 candidate genes in a 1.9 Mb genomic interval. Most of these genes carry fixed-point mutations that minimally affect their coding sequence and functional domains. Only a multi-nucleotide polymorphism in the Premelanosome protein (PMEL) results in the disruption of the protein sequence. We also characterize the different PMEL isoforms expressed in Terrazzo individuals and analyse the PMEL expression profile both with bulk RNA sequencing and real-time quantitative PCR in the developing skin and other embryonic tissues of Terrazzo and wild-type individuals. To confirm that PMEL is involved in the proper pattern establishment in corn snakes, we produce PMEL corn snake knockouts (KO), that display striped phenotypes, like the Terrazzo individuals. Electron microscopy imaging of wild-type, Terrazzo and PMEL knockouts reveal modifications in the subcellular structure of melanophores and xanthophores. Furthermore, we perform single-cell transcriptomic analyses to investigate the PMEL expression profile in chromatophores, and we show that it is expressed by all chromatophore progenitors at varying levels. Whole-mount in situ hybridisations on embryos of relevant developmental stages using a PMEL probe document how the blotches of the wild-type animals are formed. It is not possible to detect any PMEL expression on the Terrazzo embryos, either because its expression level is too low, or because the number of cells expressing it is reduced. Additional whole-mount hybridisations with different chromatophore progenitors’ markers, identified by our single-cell experiment, suggest that the overall number of chromatophores is reduced in Terrazzo individuals. Thus, a minimum number of chromatophores is necessary for the formation of blotches in the wild-type. These results will serve as a starting point to elucidate the developmental mechanisms involved in reptilian colouration and have implications in the fields of evolutionary developmental biology, herpetology and ecology.Fig. 1 Skin colour pattern of wild-type and Terrazzo corn snakes.

Dorsal overview (left), close-ups of lateral (top centre) and ventral views (bottom centre), and head dorsal view (right) of a wild-type (a) and a Terrazzo (b) individual. Close-ups of the eyes are provided in the insets.

Results

Description of the Terrazzo corn snake morph

Wild-type corn snakes are characterised by red dorsal and lateral blotches, delineated by a black line, on an orange background and the black pupil of their eyes is surrounded by a red iris (Fig. 1a). The Terrazzo corn snake morph originates from a cross of two 'rosy rat' corn snakes. These 'rosy rat' corn snakes are natives of the Key islands in Florida and they are hypomelanistic, with reduced black pigmentation17. The Terrazzo hatchlings have two red dorsal stripes that run down the entire length of their body on an orange background (Supplementary Fig. 1). Two red lateral stripes are visible at the anterior part of their body, but these gradually fade towards the posterior part. The longitudinal stripes are delineated by very thin black lines, whose thickness varies along the length of the body and from one individual to the other; occasionally, they can be absent. Over the years and as the animal grows, the dorsal and lateral stripes are maintained at the anterior part of the body, and the posterior part takes up an overall speckled grey/orange colouration without a distinct pattern (Fig. 1b). The iris of the Terrazzo animals is also speckled, instead of red. Finally, the black and white checkers on the ventral scales of the wild-type are replaced by a uniform white colouration at the anterior part of the body of the Terrazzo animals that progressively turns red posteriorly (Supplementary Fig. 1). Note that red pigmentation can also be seen ventrally at the posterior body of the wild-type animals.

Mapping of the Terrazzo variant on PMEL

Crosses in our facility confirmed that a single-locus recessive mutation is responsible for the Terrazzo morph18. We performed mapping-by-sequencing analyses to retrieve the causative variant, as previously described19. In short, we crossed a homozygous Terrazzo male (tz/tz genotype) with two heterozygous females (tz/+) to obtain homozygous and heterozygous offspring for the causative allele. We sequenced four genomic DNA libraries from: the homozygous male, the two heterozygous females, a pool of 26 homozygous offspring and a pool of 19 heterozygous offspring (Supplementary Table 1; NCBI accession number PRJNA1073516). We aligned each library separately to our corn snake genome assembly19 (GCF001185365.1) and searched for single- and multi-nucleotide polymorphisms co-segregating with the Terrazzo genotype in non-repetitive elements19. We could identify an interval of 12.05 Mb on Super-scaffold 85 (NW_023010713; from 2.86 to 14.91 Mb; Fig. 2a, b), where there are 42,973 co-segregating SNPs/MNPs with a density of 3568 variants/Mb. We found regions with a high proportion of co-segregating variants in four additional Super-scaffolds, but not as high as in Super-scaffold 85 (Supplementary Fig. 2). We genotyped 91 individuals and identified 13 recombinants that allowed us to reduce the interval to 1.9 Mb (from 13 to 14.9 Mb; Supplementary Fig. 3a, b). The corresponding interval on the latest chromosome-length corn snake assembly from Columbia University20 (CU assembly; GCF_029531705.1) is 2.12 Mb long (NW_026844023.1; from 82.93 to 84.85 Mb; Supplementary Fig. 3c).Fig. 2 Mapping of the Terrazzo variant.

a Proportion (y-axis) of quality-filtered biallelic SNP/MNP co-segregating with the Terrazzo locus in the four genomic libraries compared to informative quality-filtered parental variants. Proportions are calculated for scaffolds >1 Mb, with a 1-Mb sliding window and a step of 100 Kb. Scaffolds (alternatively coloured black and red) are ordered from longest to shortest. We label Super-scaffold 85 and the other four super-scaffolds with a high proportion of co-segregating SNPs (asterisks). A close-up of these super-scaffolds is provided in Supplementary Fig. 2. b Proportion of biallelic variants (SNP/MNP and indels) in 200-Kb intervals with a 50-Kb step co-segregating with the Terrazzo locus in Super-scaffold 85. Dark green dots correspond to the 12 Mb region with the highest proportion of cosegregating variants. The shaded green area corresponds to the reduced 1.9 Mb interval based on genotyping of additional individuals (Supplementary Fig. 3). An arrow points to the position of PMEL. c A nucleotide polymorphism at the end of exon 8 and the beginning of intron 8 affects the PMEL transcript in Terrazzo individuals (Supplementary Fig. 4). d InterProScan domains of the wild-type PMEL. The alternative Terrazzo PMEL proteins are 478 and 532 amino acids long (red vertical lines). The red arrows point to the target sites of the gRNAs used for gene-editing with CRISPR-Cas9. The range of the two fragments, Mα and Mβ, after cleavage are shown below. IP signal IP, CAF core amyloid fragment, PKD polycystic kidney disease domain, RPT repeat domain, CS proprotein convertase cleavage site, KLD Kringle-like domain, TM transmembrane helix, cyto cytoplasmic domain.

Based on the NCBI annotation of both assemblies, there are 35 protein-coding genes in the reduced interval (Supplementary Table 2). We identified indels, single- and multi-nucleotide polymorphisms that co-segregate with the Terrazzo genotype in 27 of these genes, most of which result in synonymous substitutions. There are non-synonymous amino acid modifications in the coding sequence of 14 genes, but these do not affect the domain predictions by InterProScan21. On the other hand, a multi-nucleotide polymorphism, from ACG to TCC (Super-scaffold 85; position 13,272,041) falls at the end of exon 8 (one nucleotide – A to T) and on the donor site of intron 8 of the Premelanosome Protein (PMEL) gene (two nucleotides – CG to CC) (Fig. 2c).

PMEL is a glycoprotein produced in melanophores and the retinal pigmented epithelium. The PMEL protein goes through substantial post-translational modifications, mainly O-glycosylation and cleaving, during and after transfer in the early-stage melanosomes. There, the PMEL fragments assemble into striated amyloid fibrils providing a scaffold for the deposition of melanin22. The corn snake PMEL protein is 713 amino acids (aa) long and comprises: an initial signalling peptide at the N-terminus (IP; aa 1–20), the core amyloid fragment (CAF; aa 152–227), the polycystic kidney disease domain (PKD; aa 238-299), the repeat domain (RPT; aa 306–510), the Kringle-like domain (KLD; aa 559–609), a transmembrane helix (TM; aa 638–660), and a cytoplasmic domain (cyto; aa 661–713)21,23–25 (Fig. 2d). The CAF and the RPT domains are not identified by InterProScan, but they have previously been characterised in the corn snake PMEL25. Cleavage at the proprotein convertase site (CS; aa 506–511) produces the lumenal Mα fragment, which is disulphide-linked to the membrane-integrated Mβ fragment via the KLD domain. The amyloid core is produced by further fragmentation of the Mα fragment26, with the interaction of the CAF and the PKD domains. The RPT domain is essential for proper amyloid morphology to maintain the sheet architecture25.

Two PMEL isoforms have been annotated on the corn snake genome (XM_034412587 and XM_034412588) that differ by three nucleotides; the first is 3135 bp long and the second 3132 bp. In the second isoform, the three first nucleotides of exon 8 are missing. To verify how the multi-nucleotide polymorphism identified in Terrazzo impacts the transcript, we sequenced the PMEL transcript from one homozygous and one heterozygous Terrazzo individual (Supplementary Fig. 4a) using mRNA extracted from the dorsal skin of embryos at embryonic day 20 (E20; corresponding to the developmental stage 7 described in the Supplementary text). We identified three isoforms from the homozygous individual: (i) isoform TZ1 missing exon 8, (ii) isoform TZ2 with a 20-nucleotides insertion between exons 8 and 9, which corresponds to the first 20 nucleotides of intron 8, and (iii) isoform TZ3 with a 233-nucleotides insertion between exons 8 and 9, which corresponds to the first 233 nucleotides of intron 8 (Supplementary Fig. 4b). The three isoforms produce truncated PMEL proteins—the TZ1 protein is 478 aa long and the TZ2 and TZ3 proteins are 532 aa long—that lack the Mβ fragment, which includes the Kringle-like domain and the transmembrane helix. We cannot exclude the presence of other isoforms. We could only sequence the WT isoform from the heterozygous individual, but the other isoforms are likely to be expressed at lower levels not detectable with our amplification and sequencing method.

PMEL expression is reduced in Terrazzo corn snakes

Next, we performed differential expression analyses by bulk RNA sequencing to verify the expression levels of PMEL during development. We extracted RNA from the skin of three homozygous (tz/tz) and three wild-type (+/+) embryos at E20 (developmental stage S7; Supplementary text). This is the earliest stage at which the skin can be dissected, and early-stage unpigmented chromatoblasts—progenitors of melanophores, xanthophores and iridophores—are already present in this tissue. Note that the six embryos were from the same cross (tz/+ x tz/+) and clutch to minimise expression differences due to a variable genetic background and developmental stage (Supplementary Fig. 5a, b). The most significantly downregulated gene is PMEL (adjusted p value: 4.2326 × 10−165; log2 fold-change: −3.87; Fig. 3a). Among the 12 genes with an absolute log2 fold-change greater than 1, only PMEL is within the genomic interval of the Terrazzo mutation (Fig. 3b). Transcripts per million (TPM) calculations at the exon level showed that all exons are less expressed in Terrazzo compared to the wild-type (Supplementary Fig. 5c). Thus, not only the structure, but also the expression level of PMEL is affected in the developing skin. There are three possibilities for this observation: (i) an additional undetermined mutation of a PMEL regulatory element alters its expression levels, (ii) the RNA surveillance mechanism is activated by the premature STOP codon and initiates the rapid degradation of the mutated PMEL mRNA27 by nonsense-mediated mRNA decay, and (iii) the number of cells expressing PMEL is reduced.Fig. 3 PMEL expression levels.

a Volcano plot depicting statistically significant gene expression changes between wild-type and homozygous Terrazzo embryonic dorsal skin samples in terms of log2 fold-change (x-axis) and negative log10 of p value (y-axis). In dark cyan, the genes significantly upregulated in Terrazzo and in Indian red, the ones significantly downregulated (cut-off adjusted p value: 0.05 based on a Wald test and an FDR correction). We indicate the name of the differentially expressed genes situated in the reduced genomic interval, or nearby. b Heatmap of genes with an absolute log2 fold-change greater than 1. PMEL is highlighted in red and APOH is in bold, because it is situated near the genomic interval. c Gene expression ratio of PMEL in the embryonic dorsal skin, the brain and a body part using two different primers sets, the first spanning exons 3 and 4 (left bar plot of each pair with a black border) and the second exons 8 to 10 (right bar plot of each pair without a border). The PMEL expression for the Terrazzo samples is barely detectable (2−(ΔΔCt) < 0.0003), when using the second set of primers spanning the mutated site. The number of individuals used for each category is indicated next to the corresponding bars. Data are presented as mean ± s.d. Source data are provided as a Source Data file.

ERBB3 (Receptor tyrosine-protein kinase erbB-3), a gene involved in the differentiation of metamorphic melanophores in zebrafish28, is present in the interval, but we do not detect any amino acid modifications in its coding sequence (Supplementary Table 2). It is significantly downregulated at E20 (adjusted p value <0.05), but the log2 fold-change and its expression level are very low (Fig. 3a and Supplementary Fig. 5d). Also, the Apolipoprotein H (APOH) gene is located near the interval, and it presents a similar expression profile to ERBB3 (Fig. 3a and Supplementary Fig. 5d). As these two genes are within or near the genomic interval, we cannot exclude their involvement in the Terrazzo phenotype, despite their low levels of expression both at E20/S7 and at E25/S9 (Supplementary Fig. 5d, e).

To confirm the results of the bulk RNA-seq analyses, we also performed real-time quantitative PCR on a greater number of individuals and from different embryonic tissues (Fig. 3c). We extracted RNA from the skin of 27 E20 embryos produced by three different tz/+ x tz/+ crosses. We used two sets of primers, the first spans exons 3 to 4 and the second exons 8 to 10. With the first set of primers, we observe only 5.5% PMEL expression in homozygous (tz/tz) individuals compared to the wild-type ones (+/+) and 55% in the heterozygous. With the second set of primers, practically no expression (0.006%) is observed in the homozygous samples, because the forward primer spans the 8/9 exon boundary, where the MNP is situated. To verify if the reduced levels of PMEL expression is specific to the skin, we also analysed RNA from the brain and from a body piece consisting of connective tissue, developing muscles and bones (no skin and no internal organs) of eight additional individuals (three tz/tz and five tz/+). We observe the same tendency; lower levels of PMEL expression in the brain and the body of homozygous animals compared to the wild-type (Fig. 3c). There are two possible and non-exclusive explanations for these observations: (i) as previously postulated, the presence of an additional undetermined regulatory mutation, and (ii) the reduced number of cells expressing PMEL.

CRISPR-Cas9-mediated inactivation of PMEL in snakes

Disruptive and missense mutations of PMEL have been reported in numerous animals (for example, mouse29, horse30, dog31, chicken32, and zebrafish33) and, in general, result in hypopigmentation, altered colouration patterning and defective melanosome biogenesis. Our genomic and transcriptomic analyses imply the presence of structural and possibly regulatory mutations affecting the Terrazzo PMEL. We thus proceeded to functional analyses and, using our established gene-editing protocol in corn snakes34, we generated PMEL knockouts. In short, we selected two gRNAs that successfully produced PMEL-mutated corn snake cultured fibroblast cell lines by targeting exons 2 (204f) and 3 (302r) of PMEL, respectively (Fig. 2d). The gRNAs were separately injected into the oocytes of three wild-type females (two with 204 f and one with 302r), which were crossed with wild-type males. In total, we injected 59 oocytes and the females laid 42 eggs, 30 of which hatched to give 19 mutated animals (Supplementary Table 3 and Supplementary Fig. 6). As was the case for the EDARADD mutants34, most induced mutations impacted both the paternal and the maternal alleles. Only one individual is heterozygous with a different mutation on each allele. The characterisation is performed on gDNA extracted from skin sheds, which is of epidermal origin, so we cannot exclude genetic mosaicism in other tissues.

Nine mutants have a wild-type skin colouration phenotype (Fig. 4a). The induced mutations in these individuals minimally impact the coding sequence of PMEL (e.g. insertion/deletion of one or two amino acids). The remaining ten carry disruptive mutations (e.g. premature STOP codons), and their skin phenotype is affected. Of these, seven mutants present dorsal and lateral stripes, instead of blotches, and completely lack the black checkers on the ventral side, a phenotype that strongly resembles that of Terrazzo (Fig. 4b, c). Another individual presents stripes at the anterior part and blotches at the posterior part (Fig. 4d), whereas the other two have modified blotches (Fig. 4e). These three individuals carry similar mutations as the ones with stripes. Sequencing DNA from skin biopsies taken at different regions from these animals supports non-mosaic status for this tissue (Supplementary Fig. 6). Their variable phenotype might be due to their genetic background, undetected off-target mutations, or variation in the expression levels of PMEL during their development. Note that the cells sampled with the biopsies are not the ones that established the pattern during development.Fig. 4 Gene-editing of PMEL results in stripes instead of blotches.

Dorsal, lateral and ventral views of a a corn snake with a wild-type phenotype, b a corn snake with a Terrazzo phenotype, c a PMEL mutant with a striped phenotype, d a PMEL mutant with stripes at the anterior part of the body and blotches at the posterior and e a PMEL mutant with blotches that are not properly formed. In the lateral view of (b, c), we added photos of the posterior part of the animals (on the right of the vertical dashed lines), because the pattern changes compared to the anterior part. All animals were four months old when imaged.

We observed subtle differences between the PMEL mutants and the Terrazzo animals. The dorsal and lateral stripes of the PMEL mutants are clearly delineated by black lines, whereas the black border is reduced or even absent in Terrazzo animals. This difference could be due to the hypomelanistic background of the Terrazzo individuals. Indeed, the Terrazzo mutation first appeared in ‘rosy rat’ corn snakes, a hypomelanistic morph that exists in nature. As it has already been shown in the zebrafish1, several genes are expected to be involved in the patterning of corn snakes chromatophores during development, and PMEL is only one of them. Furthermore, it remains to be seen if, over the years, the colouration of the PMEL mutants will become speckled at the posterior part, as is the case for the Terrazzo adults. It will also be necessary to check the expression levels of PMEL in the offspring of the PMEL mutants once available (generation time of 4 years).

Subcellular modifications in the chromatophores of Terrazzo and PMEL mutants

The causative mutation of the fading vision (fdv) zebrafish mutant is a premature STOP codon in PMEL. The predicted truncated protein lacks part of the repeat domain, the Kringle-like domain, the transmembrane helix, and the cytoplasmic domain33. This modification closely resembles the one predicted for the Terrazzo PMEL. The fdv zebrafish are characterised by a reduced number of melanophores, whose amount of pigmentation and subcellular morphology is affected. Indeed, electron microscopy imaging showed that the melanosomes of fdv zebrafish are smaller, compared to the zebrafish wild-type melanosomes, and at variable stages of maturation. Our transmission electron microscopy (TEM) imaging revealed that the subcellular structure of both melanophores and xanthophores is greatly impacted in Terrazzo individuals from an early age compared to the wild-type (Fig. 5a–d). Indeed, melanosomes are smaller with an irregular, instead of elliptical, shape and the xanthosomes are filled with amorphous material, rather than concentric lamellae. Note that the melanosomes take their elliptical shape only when the amyloid fibrils formed by PMEL are properly assembled. Due to the irregular shape and variable size of the guanine crystals in the wild-type, it is not possible to evaluate if iridophores are impacted as well in the Terrazzo. However, no obvious changes in their arrangement and density were observed. The subcellular structure of the melanophores is also affected in the PMEL mutants with small melanosomes of irregular shape (Fig. 5e). Some xanthosomes seem to contain concentric lamellae, but most of them only have amorphous material.Fig. 5 TEM imaging reveals subcellular differences among wild-type, Terrazzo and PMEL knockout corn snakes.

Micrographs of melanophores, xanthophores and iridophores from the dorsal skin of a a wild-type juvenile, b a wild-type adult, c a Terrazzo juvenile, d a Terrazzo adult and e a PMEL knockout (KO) juvenile. White and green dashed circles highlight normal and modified melanosomes, respectively. Black and cyan dashed circles highlight xanthosomes with concentric lamellae and amorphous material, respectively. Three blocks of each type of sample were prepared. Scale bars: 1 μm.

In the wild-type, there are melanophores and xanthophores within the blotches, and xanthophores and iridophores in the background with only scarce melanophores19. In the adult Terrazzo skin, we found scarce epidermal and dermal melanophores and numerous xanthophores within the stripes, whereas only epidermal melanophores and iridophores were present in the background (Supplementary Fig. 7a, b). Xanthophores are most likely present in the background as well, but in small numbers. Based on this arrangement of the chromatophores, we assume that the Terrazzo stripes roughly correspond to the wild-type blotches. At the posterior part of the body with speckled colouration, all types of chromatophores are present (Supplementary Fig. 7c). In the PMEL mutant, we also mainly find xanthophores in the stripe and iridophores in the background (Supplementary Fig. 8).

PMEL is expressed by all chromatophore progenitors

Based on single-cell RNA sequencing experiments in zebrafish, pmela—one of the two pmel copies in this species—is mainly expressed in melanophores and their progenitors5. To investigate if this is the case for corn snakes, we performed single-cell RNA sequencing (Fig. 6) on dissociated cells from the dorsal skin of a wild-type S9 embryo (at E25). We staged the embryos as described in the Supplementary text and based on the staging of African house snake embryos35. Although no pigmentation is visible on the S9 embryo, the development of the skin is well advanced with scale placodes formed and the patterning of the chromatophores underway (Fig. 7). After quality check, pre-processing, and filtering, we could analyse the transcriptome of 3570 cells (Supplementary Fig. 9a, b). We clustered the cells in an unsupervised manner, and we identified nine distinct cell clusters, two of which comprise 91.3% of the cells (Fig. 6a). The expression of collagen genes by these cells suggests that they belong to two fibroblast populations (Supplementary Fig. 9c). Indeed, comparison with early mouse skin markers of fibroblasts36,37 (Supplementary Fig. 9d, e) showed that a significant proportion of upregulated transcripts in cluster 1 overlap with markers of the upper dermis and dermal condensate fibroblasts. The significantly upregulated transcripts in cluster 2 correspond to markers of fascia-forming and pre-adipocyte markers, as well as to markers of muscle-supportive fibroblast. Based on the differentially expressed genes of the remaining clusters (Supplementary Data 1), we could identify them as: chromatophores (SOX10 and PMEL), epidermal cells (keratin, type I cytoskeletal 24-like gene (LOC117664631) and keratin, type II cytoskeletal 5-like gene (LOC117662708)), radial glial cells (SOX10, FABP7 and SOX2), skeletal muscle cells (MYOD1), smooth muscle cells (ACTA2), endothelial cells (TEK), and erythrocytes (haemoglobin subunits alpha-A, alpha-D, beta-1 and beta-2)38.Fig. 6 Single-cell RNA sequencing analysis of the developing corn snake skin.

a Cell-type assignment of the nine clusters produced by unsupervised clustering. b Dot plot displaying the mean expression of the top 12 genes of the chromatophore cluster in all clusters, as well as the fraction of cells expressing them in each cluster. c Gene expression levels of selected genes in the cells of the chromatophores cluster.

Fig. 7 Whole-mount in situ hybridisation of wild-type embryos with a PMEL probe.

Dorsal view of the anterior and mid-body part of embryos at developmental stages a S4, b S5, c S6, d S7, e S8 and f S9. The stages are described in the Supplementary text. At least two embryos were prepared for each stage. Scale bars: 1 mm.

Among the top 12 upregulated genes in the chromatophores cluster, we find SOX10 (SRY-Box Transcription Factor 10), which is also upregulated in glial cells, reflecting the neural crest origin of both cell types (Fig. 6b). Most of the other genes are predominantly expressed in the chromatophores cluster and they are known to be involved in melanophore development and physiology, for example, CDH3 (P-cadherin), whose expression promotes intrafollicular melanogenesis39, MLPH (Melanophilin), which participates in melanosome transportation40, and IGSF11 (Immunoglobulin Superfamily Member 11), that is required for the migration and survival of melanophores41. Only S100A1 (S100 Calcium Binding Protein A1) and TP53I11 (Tumour Protein P53 Inducible Protein 11) show a wider expression in different clusters. S100A1 was previously shown to be expressed in the brain and in skeletal and heart muscles, where it is involved in cardiac contractility regulation42. TP53I11, a target of TP53, has a tumour-suppression function43 and a ubiquitous expression44,45.

Within the chromatophores cluster (Fig. 6c), we observe a broad expression of MITF (Melanocyte Inducing Transcription Factor), which is consistent with previous findings in the zebrafish5. Indeed, in zebrafish, mitfa is expressed by all pigment progenitors and persists in subpopulations of melanophores and xanthophores5, despite that it is only required for melanophore fate specification46. Our previous findings in a closely related species, the Texas rat snake, suggest that, in snakes, MITF plays a role in the differentiation of both melanophores and xanthophores47. Indeed, MITF Texas rat snake mutants are leucistic, lacking melanophores and xanthophores, but maintaining iridophores. We thus assume that at this stage of development, the chromatophores cluster mainly consists of progenitors at early stages of differentiation. PMEL is expressed by all chromatophores at varying levels. These progenitors also express early melanophore differentiation markers GPNMB (Glycoprotein Nmb − 95.1% cells), and MLANA (Melan-A – 94.12% cells). On the other hand, late differentiation melanophore markers, such as TYR (Tyrosinase), are expressed by a small subset of cells and in low levels. Similarly, only a few cells express xanthophore and iridophore determination markers, such as (i) PAX7 (Paired Box 7), involved in zebrafish xanthophore differentiation48, (ii) SOX5 (SRY-Box Transcription Factor 5), which is implicated in xanthophore specification in zebrafish and leucophore specification in medaka49, and (iii) LTK (Leucocyte Receptor Tyrosine Kinase), expressed by zebrafish iridophores and their precursors50. Presumably these cells have committed to the xanthophore and iridophore lineages, despite still expressing melanophore markers. Note that no pigmentation is visible at the developmental stage we sampled. Melanin production only starts after E32, and the iridophores and xanthophores fully mature after E40.

Spatio-temporal expression pattern of PMEL in developing snakes

Next, we performed whole-mount in situ hybridisations (WISH) using a PMEL probe to elucidate its expression pattern in the developing skin. In the wild-type, PMEL is first expressed at S3 in individual cells that form two dorsal stripes (Supplementary Fig. 10a). At this stage, the skin over the neural tube is very thin and the coelum is open ventrally, so we consider that the left and right sides of the embryo develop independently34. At S4 and S5 (Fig. 7a, b), a greater number of cells express PMEL. At the anterior part near the head, the skin over the neural tube becomes thicker and is populated by labelled cells. At the mid-body part, the labelled cells start to form aggregates, which are not aligned on the left and right sides. At S5, scale placodes become apparent, as regions devoid of labelled cells. At S6 (Fig. 7c), the skin over the neural tube is fully developed, and the scale placodes cover the embryonic skin entirely. At this stage, the aggregates of PMEL-positive cells are not easily discernible, possibly due to the staining. From S7 to S9 (Fig. 7d–f), the aggregates of labelled cells gradually become more defined and finally resemble the dorsal blotches of the adult (Fig. 1a), with the alignment of the left and right side. Note that labelled cells can also be found in the background, around the future blotches, but their concentration is reduced. Finally, we observe that at these stages of development, the chromatophore progenitors avoid the regions of the scale placodes, in a similar manner as chromatophores avoid the developing lateral line in salamanders51. We also performed WISH on a developmental series of Terrazzo embryos, and we were not able to detect any cells expressing PMEL (Supplementary Fig. 10b, c), despite that the probe spans exons 2 to 6, an intact region of the Terrazzo PMEL transcript. The PMEL expression, if any, is likely below the detection threshold of the hybridisation staining in the Terrazzo individuals.

To elucidate the impact of the PMEL mutation on the number of chromatophores in Terrazzo, we performed WISH with an MLANA and a GPNMB probe, because they are expressed by the majority of the chromatophores at similar levels as PMEL based on our single-cell transcriptomic analyses (Fig. 6c). Thus, if the number of chromatophores remains the same in Terrazzo, we expect the expression levels of MLANA and GPNMB to be similar in wild-type and Terrazzo individuals, independently of the patterning process. We rather observe a very low staining with probes for both genes in Terrazzo (Fig. 8), at the regions where the dorsal stripes will eventually form. Although accurate quantification is not possible with WISH, in the higher magnification images, we see that the number of positive cells is lower in Terrazzo, suggesting a reduction in the number of chromatophores. Note that reduced expression of these markers in Terrazzo individuals is also confirmed by the bulk RNA-seq data (Supplementary Fig. 5f). These observations suggest an overall reduction of the chromatophore cells in Terrazzo, rather than the presence of an undetermined regulatory mutation that would exclusively affect the expression of PMEL. Furthermore, it seems like the patterning process is arrested in Terrazzo individuals and the aggregates and blotches do not form, thus the presence of stripes.Fig. 8 Whole-mount in situ hybridisation of wild-type and Terrazzo embryos with MLANA and GPNMB probes.

Dorsolateral view of the mid-body part of wild-type and Terrazzo embryos at different developmental stages on which we performed WISH with MLANA (a) and GPNMB (b) probes. The regions in the dashed rectangles are magnified at the bottom. At least two embryos were prepared for each stage. Scale bars: 1 mm; magnified areas: 0.5 mm.

Based on the single-cell transcriptomic data and the WISH experiments, we conclude that the subtle differences in the expression profile of the chromatophore progenitors seem sufficient for the early establishment of the blotched pattern (Fig. 7). It is likely that most cells expressing PMEL at a detectable level form blotches, similarly to the adult pattern. The expression of PAX7, a classical xanthophore marker, is barely detectable by whole-mount in situ hybridisations, but it seems to be mainly present in the background region (Supplementary Fig. 11a). This result should be taken with caution and highlights the need to identify species-specific xanthophore and iridophore markers, given that the markers from zebrafish studies, such as CSF1R52 and GCH248 for xanthophores, as well as TFEC53 and CSF152 for iridophores, are barely expressed in corn snakes at this stage of development when the pattern is established (Supplementary Fig. 11b).

Discussion

Our genomic, bulk, and single-cell transcriptomic, functional, and developmental characterisation of the Terrazzo corn snake colour morph and the extensive comparison with wild-type corn snakes puts forward the dual role of PMEL in snake skin colouration, both in the patterning of chromatophores during embryogenesis and the melanogenesis in melanophores. Our subcellular imaging revealed that the Terrazzo melanophores are smaller and misshaped, probably due to the truncated PMEL protein they produce. Indeed, in the absence or modification of the Mβ fragment (including the Kringle-like domain and the transmembrane helix), it is possible that PMEL is not properly transferred inside the early-stage melanosomes. The presence of PMEL amyloid fibrils formed by the Mα fragment in the cytoplasm could be toxic for the cells. This toxicity could explain the gradual transition from a striped to a speckled pattern in adult Terrazzo animals. We were intrigued to also find modifications in the subcellular structure of Terrazzo xanthosomes, which were missing concentric lamellae. Indeed, very little is known about the formation of the lamellae in xanthosomes, and PMEL has never been associated with them. These results should be taken with caution as they could be linked to the genetic background of the Terrazzo animals, rather than the Terrazzo causative mutation. It is worth noting, though, that, besides the melanosomes, the xanthosomes of the PMEL mutants we produced were also impacted. We now need to wait for these gene-edited animals to grow, such that we can verify if the subcellular structure of their chromatophores and their skin pattern will change over time.

Both the Terrazzo corn snakes and the PMEL mutants we generated are characterised by stripes rather than blotches, confirming the involvement of PMEL in the patterning process. From our single-cell experiments and developmental studies, we conclude that PMEL is a marker of chromatophore progenitors, which we can use to follow the pattern formation during development. In wild-type embryos, the chromatophore progenitors homogeneously populate the developing skin, while they differentiate from neural crest cells. As their numbers increase, they start to independently form aggregates on the left and right sides of the embryo. When the skin of the two sides merges, the aggregates align, and the future blotches become apparent. Note that we also observed such a uniform distribution of PMEL-expressing cells during the establishment of the leopard gecko skin pattern14, before the bands of the hatchlings are formed. In Terrazzo embryos, expressing low levels of mutated PMEL isoforms, we speculate that the patterning process is arrested at an early stage, when the PMEL-expressing cells are evenly distributed at each side of the dorsal midline, resulting in stripes rather than blotches. It is possible that the number of chromatophore progenitors in Terrazzo is reduced and, therefore, the aggregates do not form. Thus, we postulate that a minimum number of chromatophore progenitors is necessary for the aggregates and, later on, the blotches to form.

In Nile Tilapia fish, double pmela and pmelb knockouts have a smaller number of melanophores and a greater number of xanthophores54. In the case of the Terrazzo corn snakes, we suspect that the number of all chromatophore progenitors is reduced. Alternatively, changes in PMEL could primarily impact melanophores differentiation and consequently the patterning of the other chromatophores, in a similar manner as the overexpression of csf1a primarily impacts xanthophores, but then the patterning of all chromatophores is modified in zebrafish55. It is possible that when removing the dorsal skin, we primarily sample the melanophore progenitors, so sampling different parts of the embryo could reveal the presence of progenitors of the other lineages. Note that our single-cell transcriptomic analyses on the dorsal skin of a wild-type leopard gecko embryo also support a uniform expression of PMEL by all chromatophore precursors, before they differentiate and mature14.

To support our hypothesis on pattern formation, we would need to perform single-cell experiments at different developmental stages and describe the dynamic expression profiles of the chromatophore progenitors as they differentiate and mature. These transcriptomic experiments would also allow us to identify markers of each type of chromatophores and investigate their spatial distribution at each developmental stage. The establishment and optimisation of gene-editing experiments will also enable us to understand the role of key players in the colouration of other vertebrates, such as ERBB3, whose expression level we found low in corn snake and its function in reptilian colouration, if one, remains unknown. As we continue the characterisation of different colour and colour pattern morphs of corn snakes and other squamate model species, we will deepen our understanding of the developmental processes that bring about the spectacular diversity of their skin colouration and unravel the evolution of the molecular mechanisms responsible for animal colouration.

Methods

Experimentation model

Corn snakes were housed and bred at the LANE animal facility running under veterinary cantonal permit no. 1008. Sampling and imaging were performed under the experimentation permits GE24/33145 and GE150/34215.

Crosses and deep sequencing of Terrazzo individuals

Genomic DNA was extracted from the parents and the offspring using the QIAGEN DNeasy Blood and Tissue kit (69504) following the manufacturer's instructions. The TruSeq DNA PCR Free libraries were sequenced using an Illumina HiSeqX instrument, producing 151 bp paired-end reads (Macrogen). We obtained 173.1 to 185.9 million paired-end reads per library. We checked the data quality and the presence of adaptors with FASTQC. We performed a quality filtering with sickle v1.3356. We retained between 167.4 and 179.8 million reads, which correspond to a 29.7-31.9x average coverage for a 1.7 Gb genome.

Variant calling

The four genomic libraries (two for the parents and two for the offspring) were aligned to the corn snake genome19 using bwa v0.7.1657 with default parameters in mem mode. We used SAMtools v1.958 to: convert the output SAM files into BAM, to remove duplicates using the fixmate mode with the -m flag and the markdup mode with the -r flag, and to sort out the reads by their leftmost coordinates. We identified genomic variants with Platypus v0.8.159 and retrieved the genomic interval where the Terrazzo locus is located as previously published19.

PMEL transcript amplification and cloning

We extracted mRNA from embryonic tissue with the Direct-zol RNA Miniprep kit (Zymo Research, R2081) and prepared complementary DNA with PrimeScript™ Reverse Transcriptase (Takara Bio, RR047A). We then amplified and sequenced the PMEL transcript in a Terrazzo and a heterozygous individual. The primers are provided in the Supplementary Table 4.

We extracted from an 1% agarose gel the three visible bands amplified with the primers EG_PMELc_1321-F/EG_PMELc_2090-R from a homozygous individual using the QIAquick Gel Extraction Kit (Qiagen, 28704). The extracted bands were Sanger sequenced with the same primers.

Bulk RNA-seq sampling and analysis

We first extracted genomic DNA from a piece of the tail from each E20 embryo, and we genotyped them using the primers EG_sc85-2900655-F/EG_sc85-2901018-R and EG_sc85-14891026-F/EG_sc85-14891444-R. We then extracted total RNA from the skin of three homozygous Terrazzo and three wild-type embryos using the Direct-zol RNA MiniPrep (Zymo Research, R2050). We sequenced between 18.8 and 21 million reads from each TruSeq Stranded mRNA Library. For the skin microdissection, we carefully detached the skin along the opening of the ventral side before pulling it off. Sufficient RNA was obtained from each embryo, so the samples were not pooled. All samples had an RNA integrity number (RIN) ≥9.6.

The bulk RNA-seq samples were aligned to the corn snake genome (GCF_001185365.1) with STAR (v2.7.0d) using default parameters for paired-end libraries. Gene expression quantification was performed using the featureCounts function implemented in the R package Rsubread (v2.2.1), counting uniquely mapped paired-end reads. The filterbyExpr function implemented in the edgeR package (v3.40.2) was used to filter out genes with a low number of counts for downstream analyses. Data normalisation, transformation (considering variance stabilising transformation), principal component analysis and differential expression analyses with the Wald test and an FDR of 0.05 were performed with the DESeq2 package (v1.38.3). To obtain the read count that maps to each PMEL exon with featureCounts, we generated a GTF file that included only PMEL with a unique ID per exon. We disregarded multimapping reads but considered overlapping reads that span more than one exon. For the TPM calculation, we corrected for the exon length and the library size.

Real-time quantitative PCR

We extracted RNA with the Direct-zol™ RNA MiniPrep (Zymo Research, R2050). Reverse transcription starting with 500 ng of RNA was performed with the PrimeScript RT Reagent Kit with gDNA Eraser (Takara Bio) which includes a step for the elimination of any gDNA contamination. We used the PowerUp SYBR Green Master Mix (Thermo Fisher Scientific) for the quantitative RT-PCR run on a QuantStudio 5 Real-Time PCR System (Thermo Fisher Scientific). All samples were run in triplicates, and the primers are provided in Supplementary Table 4. We verified the PCR efficiency60 for the amplification of PMEL and ALAS1, the reference gene, and performed the statistical analyses as previously described61.

Gene-editing with CRISPR-Cas9

General anaesthesia for the surgical procedure was induced by subcutaneous injection of alfaxalone (15 mg/kg) and medetomidine (0.1 mg/kg). The animals were then intubated and maintained on intermittent positive pressure ventilation with 3–7% sevoflurane in a mixture of oxygen and air. Reversal of the anaesthesia was accomplished by subcutaneous injection of atipamezole (0.5 mg/kg). Analgesia was accomplished in a multimodal approach: (i) pre-emptive analgesia was provided by subcutaneous injection of meloxicam (0.2 mg/kg) at the induction of anaesthesia, (ii) intraoperative analgesia was provided by the medetomidine and local infiltration of the skin incision site with lidocaine (4 mg/kg), and (iii) post-operative analgesia was provided by one daily subcutaneous injection with meloxicam (0.2 mg/kg) for 3 days. Antibiotic coverage was applied prophylactically (Ceftiofur CFA, 15 mg/kg, single subcutaneous injection). The surgical procedure to access the ovaries and inject the follicles was identical to the previously described34. Genomic DNA was extracted from the skin sheds of the offspring with the DNeasy Blood and Tissue Kit (69504, QIAGEN). We used the FastStart polymerase (12032902001, Sigma) to amplify the target regions. The gRNAs and the primers for the amplification of the target regions are detailed in Supplementary Table 4.

Single-cell RNA-seq analysis

E25 corn snake skin was isolated from a single embryo and dissociated in a single cell suspension according to a protocol used for chicken embryos62. Briefly, the tissue was treated with a solution constituted of 0.25% Dispase II (Sigma, D4693) and 0.25% Trypsin in Hank’s balanced solution (Sigma, M4780) at 37 °C for 18 min. Hank’s balanced solution supplemented with 0.25% BSA (Sigma, A3294) was used for Dispase II and trypsin inactivation, cell suspension washing and resuspension. Subsequently, the single-cell RNA library was prepared using the 10X Genomics Chromium Single Cell 3’ V3 kit. The library was sequenced on a HiSeq X Ten. Cellranger (v7.0.0) identified 5,431 cells. The reads were mapped against our corn snake genome assembly (GCF001185365.1) and the corn snake mitochondrial genome63 (AM236349.1). Further details about quality metrics, filtering thresholds and downstream analysis are provided in Supplementary Fig. 8b. The potential doublets were detected with Scrublet v0.2.364. The differentially expressed genes in each cluster were identified using the scanpy rank_genes_groups method with default parameters. Expression values were normalised to a uniform depth of 10,000 read counts per cell using the scanpy normalize_per_cell function and subsequently underwent a logarithmic transformation using the log1p function.

Whole-mount in situ hybridisations

We designed species-specific digoxigenin-labelled antisense riboprobes for PMEL, MLANA, GPNMB, and PAX7 (Supplementary Table 4). Embryos at different developmental stages were fixed in 4% paraformaldehyde and dehydrated in methanol. Whole-mount in situ hybridisations (WISH) were performed as previously described65. Embryos were imaged using the VHX-6000 (Keyence).

Transmission electron microscopy

Skin pieces of 1 mm2 were fixed and sectioned as previously described19. Micrographs were taken with a transmission electron microscope Philips CM100 (Thermo Fisher Scientific) at an acceleration voltage of 80 kV with a TVIPS TemCam-F416 digital camera (TVIPS GmbH). Large montage alignments were performed using the Blendmont command-line from the IMOD software66. The display levels of the images were adjusted with ImageJ to facilitate comparisons.

Reporting summary

Further information on research design is available in the Nature Portfolio Reporting Summary linked to this article.

Supplementary information

Supplementary Information

Peer Review File

Description of Additional Supplementary Files

Supplementary Data 1

Reporting Summary

Source data

Source Data

Supplementary information

The online version contains supplementary material available at 10.1038/s41467-024-51927-0.

Acknowledgements

We would like to thank Ingrid Burgelin, Carine Langrez, Adrien Debry and Florent Montange for their technical assistance. Jean Daraspe, Antonio Mucciolo and Rodrigue Peraldi participated in the TEM sample processing and imaging, which was performed at the Electron Microscopy Facility, University of Lausanne (Switzerland). The computations were run on the Baobab cluster of the University of Geneva (Switzerland). AUA received an iGE3 PhD award and the Sara Borrell grant (CD22/00027) from the Instituto Carlos III and NextGenerationEU. Funding was provided to A.C.T. by the SNF (grant 310030_204466), the HFSP (RGP0037/2022), the Ernst and Lucie Schmidheiny Foundation (10_2023), the Fonds Général de l’Université de Genève (23_28) and the Emile Plantamour Fund (2024/12).

Author contributions

Conceptualisation, funding acquisition, methodology, project administration, supervision, validation and writing—original draft: A.C.T. Data curation and visualisation: A.C.T., A.U.-A. and P.-Y.H. Investigation and writing—review and editing: A.C.T., A.U.-A., P.-Y.H. and M.K. Software: A.U.-A. and P.-Y.H.

Peer review

Peer review information

Nature Communications thanks the anonymous reviewers for their contribution to the peer review of this work. A peer review file is available.

Data availability

The sequencing data generated in this study have been deposited in NCBI under accession code PRJNA1073516 for the genomic sequencing, GSE262160 for the single-cell and GSE262159 for the bulk RNA sequencing data. Source data are provided with this paper.

Competing interests

The authors declare no competing interests.

Publisher’s note Springer Nature remains neutral with regard to jurisdictional claims in published maps and institutional affiliations.
==== Refs
References

1. Patterson LB Parichy DM Zebrafish pigment pattern formation: insights into the development and evolution of adult form Annu. Rev. Genet. 2019 53 505 530 10.1146/annurev-genet-112618-043741 31509458
Patterson, L. B. & Parichy, D. M. Zebrafish pigment pattern formation: insights into the development and evolution of adult form. Annu. Rev. Genet. 53, 505–530 (2019).31509458 10.1146/annurev-genet-112618-043741
2. Kondo, S., Watanabe, M. & Miyazawa, S. Studies of turing pattern formation in zebrafish skin. Philos. Trans. A Math. Phys. Eng. Sci.10.1098/rsta.2020.0274 (2021).
3. Milinkovitch MC Jahanbakhsh E Zakany S The unreasonable effectiveness of reaction diffusion in vertebrate skin color patterning Annu. Rev. Cell Dev. Biol. 2023 39 145 174 10.1146/annurev-cellbio-120319-024414 37843926
Milinkovitch, M. C., Jahanbakhsh, E. & Zakany, S. The unreasonable effectiveness of reaction diffusion in vertebrate skin color patterning. Annu. Rev. Cell Dev. Biol. 39, 145–174 (2023).37843926 10.1146/annurev-cellbio-120319-024414
4. Frohnhofer HG Krauss J Maischein HM Nusslein-Volhard C Iridophores and their interactions with other chromatophores are required for stripe formation in zebrafish Development 2013 140 2997 3007 10.1242/dev.096719 23821036
Frohnhofer, H. G., Krauss, J., Maischein, H. M. & Nusslein-Volhard, C. Iridophores and their interactions with other chromatophores are required for stripe formation in zebrafish. Development 140, 2997–3007 (2013).23821036 10.1242/dev.096719
5. Saunders, L. M. et al. Thyroid hormone regulates distinct paths to maturation in pigment cell lineages. Elife10.7554/eLife.45181 (2019).
6. Howard, A. G. A. et al. An atlas of neural crest lineages along the posterior developing zebrafish at single-cell resolution. Elife10.7554/eLife.60005 (2021).
7. McCluskey, B. M., Liang, Y. P., Lewis, V. M., Patterson, L. B. & Parichy, D. M. Pigment pattern morphospace of Danio fishes: evolutionary diversification and mutational effects. Biol. Open10.1242/bio.058814 (2021).
8. Roux, N. et al. The multi-level regulation of clownfish metamorphosis by thyroid hormones. Cell Rep.10.1016/j.celrep.2023.112661 (2023).
9. Tzika, A. C. On the role of TFEC in reptilian coloration. Front. Cell Dev. Biol.10.3389/fcell.2024.1358828 (2024).
10. Uetz, P. The reptile database. http://www.reptile-database.org/ (2024).
11. Saenko, S. V. et al. Amelanism in the corn snake is associated with the insertion of an LTR-retrotransposon in the OCA2 gene. Sci. Rep.10.1038/srep17118 (2015).
12. Ullate-Agote A Milinkovitch MC Tzika AC The genome sequence of the corn snake (Pantherophis guttatus), a valuable resource for EvoDevo studies in squamates Int. J. Dev. Biol. 2014 58 881 888 10.1387/ijdb.150060at 26154328
Ullate-Agote, A., Milinkovitch, M. C. & Tzika, A. C. The genome sequence of the corn snake (Pantherophis guttatus), a valuable resource for EvoDevo studies in squamates. Int. J. Dev. Biol. 58, 881–888 (2014).26154328 10.1387/ijdb.150060at
13. Guo, L. H. et al. Genetics of white color and iridophoroma in “Lemon Frost” leopard geckos. PLos Genet.10.1371/journal.pgen.1009580 (2021).
14. Ullate-Agote, A. & Tzika, A. C. The dynamic behaviour of chromatophores marks the transition from bands to spots in leopard geckos. Proc. Natl Acad. Sci. USA. 10.1073/pnas.2400486121 (2024).
15. Garcia-Elfring A Piebaldism and chromatophore development in reptiles are linked to the tfec gene Curr. Biol. 2023 33 755 10.1016/j.cub.2023.01.004 36702128
Garcia-Elfring, A. et al. Piebaldism and chromatophore development in reptiles are linked to the tfec gene. Curr. Biol. 33, 755 (2023).36702128 10.1016/j.cub.2023.01.004
16. Dao, U. M. et al. Stripes and loss of color in ball pythons (Python regius) are associated with variants affecting endothelin signaling. G3. 10.1093/g3journal/jkad063 (2023).
17. Love, K. & Love, B. Corn Snakes: The Comprehensive Owner’s Guide (Companion House Books, 2006).
18. Pritzel, C. Cornsnake Morph Guide (C. Pritzel, 2011).
19. Ullate-Agote A Genome mapping of a LYST mutation in corn snakes indicates that vertebrate chromatophore vesicles are lysosome-related organelles Proc. Natl Acad. Sci. USA 2020 117 26307 26317 10.1073/pnas.2003724117 33020272
Ullate-Agote, A. et al. Genome mapping of a LYST mutation in corn snakes indicates that vertebrate chromatophore vesicles are lysosome-related organelles. Proc. Natl Acad. Sci. USA 117, 26307–26317 (2020).33020272 10.1073/pnas.2003724117
20. Hoge C Patterns of recombination in snakes reveal a tug-of-war between PRDM9 and promoter-like features Science 2024 383 eadj7026 10.1126/science.adj7026 38386752
Hoge, C. et al. Patterns of recombination in snakes reveal a tug-of-war between PRDM9 and promoter-like features. Science 383, eadj7026 (2024).38386752 10.1126/science.adj7026
21. Jones P InterProScan 5: genome-scale protein function classification Bioinformatics 2014 30 1236 1240 10.1093/bioinformatics/btu031 24451626
Jones, P. et al. InterProScan 5: genome-scale protein function classification. Bioinformatics 30, 1236–1240 (2014).24451626 10.1093/bioinformatics/btu031
22. Bissig, C., Rochin, L. & van Niel, G. PMEL amyloid fibril formation: the bright steps of pigmentation. Int. J. Mol Sci.10.3390/ijms17091438 (2016).
23. Berson JF Proprotein convertase cleavage liberates a fibrillogenic fragment of a resident glycoprotein to initiate melanosome biogenesis J. Cell Biol. 2003 161 521 533 10.1083/jcb.200302072 12732614
Berson, J. F. et al. Proprotein convertase cleavage liberates a fibrillogenic fragment of a resident glycoprotein to initiate melanosome biogenesis. J. Cell Biol. 161, 521–533 (2003).12732614 10.1083/jcb.200302072
24. Chrystal, P. W. et al. Functional domains and evolutionary history of the PMEL and GPNMB family proteins. Molecules. 10.3390/molecules26123529 (2021).
25. Graham, M., Tzika, A. C., Mitchell, S. M., Liu, X. R. & Leonhardt, R. M. Repeat domain-associated O-glycans govern PMEL fibrillar sheet architecture. Sci. Rep.10.1038/s41598-019-42571-6 (2019).
26. Hee JS Mitchell SM Liu X Leonhardt RM Melanosomal formation of PMEL core amyloid is driven by aromatic residues Sci. Rep. 2017 7 44064 10.1038/srep44064 28272432
Hee, J. S., Mitchell, S. M., Liu, X. & Leonhardt, R. M. Melanosomal formation of PMEL core amyloid is driven by aromatic residues. Sci. Rep. 7, 44064 (2017).28272432 10.1038/srep44064
27. Hsu SN Hertel KJ Spliceosomes walk the line Splicing errors and their impact on cellular function RNA Biol. 2009 6 526 530 10.4161/rna.6.5.9860 19829058
Hsu, S. N. & Hertel, K. J. Spliceosomes walk the line Splicing errors and their impact on cellular function. RNA Biol. 6, 526–530 (2009).19829058 10.4161/rna.6.5.9860
28. Budi EH Patterson LB Parichy DM Embryonic requirements for ErbB signaling in neural crest development and adult pigment pattern formation Development 2008 135 2603 2614 10.1242/dev.019299 18508863
Budi, E. H., Patterson, L. B. & Parichy, D. M. Embryonic requirements for ErbB signaling in neural crest development and adult pigment pattern formation. Development 135, 2603–2614 (2008).18508863 10.1242/dev.019299
29. Martínez-Esparza MM The mouse silver locus encodes a single transcript truncated by the silver mutation Mamm. Genome 1999 10 1168 1171 10.1007/s003359901184 10594241
Martínez-Esparza, M. M. et al. The mouse silver locus encodes a single transcript truncated by the silver mutation. Mamm. Genome 10, 1168–1171 (1999).10594241 10.1007/s003359901184
30. Brunberg, E. et al. A missense mutation in PMEL17 is associated with the Silver coat color in the horse. BMC Genet.10.1186/1471-2156-7-46 (2006).
31. Clark LA Wahl JM Rees CA Murphy KE Retrotransposon insertion in SILV is responsible for merle patterning of the domestic dog Proc. Natl Acad. Sci. USA 2006 103 1376 1381 10.1073/pnas.0506940103 16407134
Clark, L. A., Wahl, J. M., Rees, C. A. & Murphy, K. E. Retrotransposon insertion in SILV is responsible for merle patterning of the domestic dog. Proc. Natl Acad. Sci. USA 103, 1376–1381 (2006).16407134 10.1073/pnas.0506940103
32. Kerje S The Dominant white, Dun and Smoky color variants in chicken are associated with insertion/deletion polymorphisms in the PMEL17 gene Genetics 2004 168 1507 1518 10.1534/genetics.104.027995 15579702
Kerje, S. et al. The Dominant white, Dun and Smoky color variants in chicken are associated with insertion/deletion polymorphisms in the PMEL17 gene. Genetics 168, 1507–1518 (2004).15579702 10.1534/genetics.104.027995
33. Schonthaler HB A mutation in the silver gene leads to defects in melanosome biogenesis and alterations in the visual system in the zebrafish mutant Dev. Biol. 2005 284 421 436 10.1016/j.ydbio.2005.06.001 16024012
Schonthaler, H. B. et al. A mutation in the silver gene leads to defects in melanosome biogenesis and alterations in the visual system in the zebrafish mutant. Dev. Biol. 284, 421–436 (2005).16024012 10.1016/j.ydbio.2005.06.001
34. Tzika AC Ullate-Agote A Zakany S Kummrow M Milinkovitch MC Somitic positional information guides self-organized patterning of snake scales Sci. Adv. 2023 9 eadf8834 10.1126/sciadv.adf8834 37315141
Tzika, A. C., Ullate-Agote, A., Zakany, S., Kummrow, M. & Milinkovitch, M. C. Somitic positional information guides self-organized patterning of snake scales. Sci. Adv. 9, eadf8834 (2023).37315141 10.1126/sciadv.adf8834
35. Boback SM Dichter EK Mistry HL A developmental staging series for the African house snake, Boaedon (Lamprophis) fuliginosus Zoology 2012 115 38 46 10.1016/j.zool.2011.09.001 22206643
Boback, S. M., Dichter, E. K. & Mistry, H. L. A developmental staging series for the African house snake, Boaedon (Lamprophis) fuliginosus. Zoology 115, 38–46 (2012).22206643 10.1016/j.zool.2011.09.001
36. Jacob T Molecular and spatial landmarks of early mouse skin development Dev. Cell 2023 58 2140 10.1016/j.devcel.2023.07.015 37591247
Jacob, T. et al. Molecular and spatial landmarks of early mouse skin development. Dev. Cell 58, 2140 (2023).37591247 10.1016/j.devcel.2023.07.015
37. Metikala, S., Chetty, S. C. & Sumanas, S. Single-cell transcriptome analysis of the zebrafish embryonic trunk. PLos ONE. 10.1371/journal.pone.0254024 (2021).
38. Karlsson, M. et al. A single-cell type transcriptomics map of human tissues. Sci. Adv.10.1126/sciadv.abh2169 (2021).
39. Samuelov L Topobiology of human pigmentation: P-cadherin selectively stimulates hair follicle melanogenesis J. Invest. Dermatol. 2013 133 1591 1600 10.1038/jid.2013.18 23334344
Samuelov, L. et al. Topobiology of human pigmentation: P-cadherin selectively stimulates hair follicle melanogenesis. J. Invest. Dermatol. 133, 1591–1600 (2013).23334344 10.1038/jid.2013.18
40. Myung, C. H., Lee, J. E., Jo, C. S., Park, J. I. & Hwang, J. S. Regulation of melanophilin (Mlph) gene expression by the glucocorticoid receptor (GR). Sci. Rep.10.1038/s41598-021-96276-w (2021).
41. Eom, D. S. et al. Melanophore migration and survival during zebrafish adult pigment stripe development require the immunoglobulin superfamily adhesion molecule Igsf11. PLos Genet.10.1371/journal.pgen.1002899 (2012).
42. Ritterhoff J Most P Targeting S100A1 in heart failure Gene Ther. 2012 19 613 621 10.1038/gt.2012.8 22336719
Ritterhoff, J. & Most, P. Targeting S100A1 in heart failure. Gene Ther. 19, 613–621 (2012).22336719 10.1038/gt.2012.8
43. Xiao T TP53I11 suppresses epithelial-mesenchymal transition and metastasis of breast cancer cells BMB Rep. 2019 52 379 384 10.5483/BMBRep.2019.52.6.173 30940320
Xiao, T. et al. TP53I11 suppresses epithelial-mesenchymal transition and metastasis of breast cancer cells. BMB Rep. 52, 379–384 (2019).30940320 10.5483/BMBRep.2019.52.6.173
44. Uhlen M Towards a knowledge-based Human Protein Atlas Nat. Biotechnol. 2010 28 1248 1250 10.1038/nbt1210-1248 21139605
Uhlen, M. et al. Towards a knowledge-based Human Protein Atlas. Nat. Biotechnol. 28, 1248–1250 (2010).21139605 10.1038/nbt1210-1248
45. Uhlen, M. et al. Human protein atlas version 23. https://www.proteinatlas.org/ENSG00000175274-TP53I11/ (2024).
46. Lister JA Robertson CP Lepage T Johnson SL Raible DW nacre encodes a zebrafish microphthalmia-related protein that regulates neural-crest-derived pigment cell fate Development 1999 126 3757 3767 10.1242/dev.126.17.3757 10433906
Lister, J. A., Robertson, C. P., Lepage, T., Johnson, S. L. & Raible, D. W. nacre encodes a zebrafish microphthalmia-related protein that regulates neural-crest-derived pigment cell fate. Development 126, 3757–3767 (1999).10433906 10.1242/dev.126.17.3757
47. Ullate-Agote, A. & Tzika, A. C. Characterization of the leucistic Texas rat snake Pantherophis obsoletus. Front. Ecol. Evol.10.3389/fevo.2021.583136 (2021).
48. Nord H Dennhag N Muck J von Hofsten J Pax7 is required for establishment of the xanthophore lineage in zebrafish embryos Mol. Biol. Cell 2016 27 1853 1862 10.1091/mbc.e15-12-0821 27053658
Nord, H., Dennhag, N., Muck, J. & von Hofsten, J. Pax7 is required for establishment of the xanthophore lineage in zebrafish embryos. Mol. Biol. Cell 27, 1853–1862 (2016).27053658 10.1091/mbc.e15-12-0821
49. Nagao, Y. et al. Distinct interactions of Sox5 and Sox10 in fate specification of pigment cells in medaka and zebrafish. PLos Genet.10.1371/journal.pgen.1007260 (2018).
50. Lopes, S. S. et al. Leukocyte tyrosine kinase functions in pigment cell development. PLos Genet.10.1371/journal.pgen.1000026 (2008).
51. Parichy, D. M. When neural crest and placodes collide: interactions between melanophores and the lateral lines that generate stripes in the salamander Ambystoma tigrinum tigrinum (Ambystomatidae). Dev. Biol. 10.1006/dbio.1996.0115 (1996).
52. Patterson LB Parichy DM Interactions with iridophores and the tissue environment required for patterning melanophores and xanthophores during zebrafish adult pigment stripe formation PLos Genet. 2013 9 e1003561 10.1371/journal.pgen.1003561 23737760
Patterson, L. B. & Parichy, D. M. Interactions with iridophores and the tissue environment required for patterning melanophores and xanthophores during zebrafish adult pigment stripe formation. PLos Genet. 9, e1003561 (2013).23737760 10.1371/journal.pgen.1003561
53. Petratou, K., Spencer, S. A., Kelsh, R. N. & Lister, J. A. The MITF paralog tfec is required in neural crest development for fate specification of the iridophore lineage from a multipotent pigment cell progenitor. PLos ONE. 10.1371/journal.pone.0244794 (2021).
54. Wang CX Xu J Kocher TD Li MH Wang DS CRISPR knockouts of pmela and pmelb engineered a golden tilapia by regulating relative pigment cell abundance J. Hered. 2022 113 398 413 10.1093/jhered/esac018 35385582
Wang, C. X., Xu, J., Kocher, T. D., Li, M. H. & Wang, D. S. CRISPR knockouts of pmela and pmelb engineered a golden tilapia by regulating relative pigment cell abundance. J. Hered. 113, 398–413 (2022).35385582 10.1093/jhered/esac018
55. Patterson LB Bain EJ Parichy DM Pigment cell interactions and differential xanthophore recruitment underlying zebrafish stripe reiteration and Danio pattern evolution Nat. Commun. 2014 5 5299 10.1038/ncomms6299 25374113
Patterson, L. B., Bain, E. J. & Parichy, D. M. Pigment cell interactions and differential xanthophore recruitment underlying zebrafish stripe reiteration and Danio pattern evolution. Nat. Commun. 5, 5299 (2014).25374113 10.1038/ncomms6299
56. Joshi, N. A. & Fass, J. N. Sickle: a sliding-window, adaptive, quality-based trimming tool for FastQ files (Version 1.33) (2011).
57. Li H Durbin R Fast and accurate short read alignment with Burrows-Wheeler transform Bioinformatics 2009 25 1754 1760 10.1093/bioinformatics/btp324 19451168
Li, H. & Durbin, R. Fast and accurate short read alignment with Burrows-Wheeler transform. Bioinformatics 25, 1754–1760 (2009).19451168 10.1093/bioinformatics/btp324
58. Danecek, P. et al. Twelve years of SAMtools and BCFtools. Gigascience. 10.1093/gigascience/giab008 (2021).
59. Rimmer A Integrating mapping-, assembly- and haplotype-based approaches for calling variants in clinical sequencing applications Nat. Genet. 2014 46 912 918 10.1038/ng.3036 25017105
Rimmer, A. et al. Integrating mapping-, assembly- and haplotype-based approaches for calling variants in clinical sequencing applications. Nat. Genet. 46, 912–918 (2014).25017105 10.1038/ng.3036
60. Pfaffl, M. W. A new mathematical model for relative quantification in real-time RT-PCR. Nucleic Acids Res.10.1093/nar/29.9.e45 (2001).
61. Livak KJ Schmittgen TD Analysis of relative gene expression data using real-time quantitative PCR and the 2(-Delta Delta C(T)) method Methods 2001 25 402 408 10.1006/meth.2001.1262 11846609
Livak, K. J. & Schmittgen, T. D. Analysis of relative gene expression data using real-time quantitative PCR and the 2(-Delta Delta C(T)) method. Methods 25, 402–408 (2001).11846609 10.1006/meth.2001.1262
62. Williams, R. M. & Sauka-Spengler, T. Dissociation of chick embryonic tissue for FACS and preparation of isolated cells for genome-wide downstream assays. STAR Protoc.10.1016/j.xpro.2021.100414 (2021).
63. Douglas DA Janke A Arnason U A mitogenomic study on the phylogenetic position of snakes Zool. Scr. 2006 35 545 558 10.1111/j.1463-6409.2006.00257.x
Douglas, D. A., Janke, A. & Arnason, U. A mitogenomic study on the phylogenetic position of snakes. Zool. Scr. 35, 545–558 (2006).10.1111/j.1463-6409.2006.00257.x
64. Wolock SL Lopez R Klein AM Scrublet: computational identification of cell doublets in single-cell transcriptomic data Cell Syst. 2019 8 281 291.e289 10.1016/j.cels.2018.11.005 30954476
Wolock, S. L., Lopez, R. & Klein, A. M. Scrublet: computational identification of cell doublets in single-cell transcriptomic data. Cell Syst. 8, 281–291.e289 (2019).30954476 10.1016/j.cels.2018.11.005
65. Riddle RD Johnson RL Laufer E Tabin C Sonic-hedgehog mediates the polarizing activity of the ZPA Cell 1993 75 1401 1416 10.1016/0092-8674(93)90626-2 8269518
Riddle, R. D., Johnson, R. L., Laufer, E. & Tabin, C. Sonic-hedgehog mediates the polarizing activity of the ZPA. Cell 75, 1401–1416 (1993).8269518 10.1016/0092-8674(93)90626-2
66. Kremer JR Mastronarde DN McIntosh JR Computer visualization of three-dimensional image data using IMOD J. Struct. Biol. 1996 116 71 76 10.1006/jsbi.1996.0013 8742726
Kremer, J. R., Mastronarde, D. N. & McIntosh, J. R. Computer visualization of three-dimensional image data using IMOD. J. Struct. Biol. 116, 71–76 (1996).8742726 10.1006/jsbi.1996.0013
