PMEL is involved in snake colour pattern transition from blotches to stripes

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.


Supplementary text Staging of the corn snake embryos
The embryos used for transcriptomic analyses and whole-mount in situ hybridisations were staged as described here and roughly based on the staging of Boaedon (Lamprophis) fuliginosus 1 .

Stage S2
The body forms four to five tight coils.The eye is lightly pigmented.The remnants of two pharyngeal slits are visible and the frontonasal prominence has a squared shape.There are still no visible endolymphatic ducts.The tip of the lower jaw is anterior to the eye.

Stage S3
The body forms four to five tight coils up to the cloaca.Bulges are observed ventrally, where the ventral scales will form, and the coelum is open along the entire body.The heart and liver are external.The endolymphatic duct spots are visible as white dots, and the diencephalon is smaller than the mesencephalon.The lower jaw is smaller than 1/3 of the upper jaw.The eye is pigmented with the choroid fissure evident.One pharyngeal slit is visible.

Stage S4
The diencephalon has a similar size as the mesencephalon.The lower jaw has similar length as the upper jaw.The choroid fissure is not visible.

Stage S5
The eye lid/scale starts to close and the skin and bone over the brain are more opaque.The lateral placodes start to form.

Stage S6
The body forms 3.5 coils and the coelum is almost closed from the head to the heart region.The heart and liver are partially internalised.The lateral placodes are fully formed.

Stage S7
The body forms 3 coils and the coelum is closed from the head to the heart regions.The heart and liver are internalised but still visible.The dorsal placodes start to form.

Stage S8
The body is fully scaled laterally.The coelum is closed along two thirds of the body length.The heart and liver are no longer visible.The ventral scales near the head are fused.The mouth starts to open, and the tongue is visible.

Stage S9
Dorsal scales start to form near the head, and the ventral scales are fused along two thirds of the body length.The ventral scales start to become asymmetrical.The coelum is open near the cloaca and the intestine is visible.

Stage S10
The dorsal scales are formed along two thirds of the body length.All ventral scales are fused.

Supplementary Figure 2 .
Additional scaffolds presenting a high proportion of cosegregating variants.Proportion of biallelic variants (SNP/MNP and indels) cosegregating with the Terrazzo locus considering a 200-kb sliding window and a step of 50 kb in Superscaffold 65 (a), Super-scaffold 41126 (b), Super-scaffold 110 (c) and Super-scaffold 82133 (d).In none of these scaffolds do we observe regions fully cosegregating with the Terrazzo locus, as is the case for Super-scaffold 85. Supplementary Figure 3. Reduction of the genomic interval harbouring the Terrazzo variant.(a) Schematic representation of Super-scaffold 85 and the positions (in Mb) where recombination events were detected.The numbers below the scaffold detail the number of recombinants at each position.In dark green, we highlight the reduced interval.(b) List of the 13 recombinant individuals (among 91) and the genotyping results.Recombined sites are in orange.(c) Position of the reduced interval on the CU assembly.Supplementary Figure 4. Amplification of the PMEL transcript.(a) Amplification of the PMEL transcript from a homozygous and a heterozygous individual with five overlapping fragments.The DNA ladder ranges from 200 to 10,000 bp (SmartLadder SF, Eurogentec).NC: PCR negative control (b) Schematic representation of the WT PMEL transcript and the three PMEL isoforms detected in Terrazzo individuals (TZ1, TZ2, and TZ3).The red lines indicate premature STOP codons.We cannot exclude the presence of additional isoforms.Supplementary Figure 5. RNA-seq gene expression analyses.(a) Principal component analyses of the three wild-type and the three Terrazzo samples used for RNA-seq.(b) Heatmap of the 60 most differentially expressed genes located on the Z chromosome (based on Crotalus viridis Z chromosome; NCBI accession CM012323.1).Given that corn snake males are homogametic (ZZ), it is clear from this comparison that sample WT1 is the only male sampled.This is the main biological difference responsible for the clustering of the samples in the PCA presented in (a).We do not expect clustering based on the coloration at this stage, because the chromatophores are at a very early stage of differentiation and no coloration is visible.(c) Transcripts per million (TPM) calculations for the 12 exons of PMEL, corrected for exon and library size.(d) TPMs for the genes in the Terrazzo interval with an adjusted p-value <0.05.We include the APOH gene which is located near the interval, and it is significantly downregulated in Terrazzo.The gene LOC117662899 has been replaced by LOC117662900 in NCBI and it corresponds to an uncharacterized non-coding RNA.(e) Comparison of the TPMs for PMEL, ERBB3, and APOH as calculated from the bulk RNAseq at E20 (S7) and the CPMs as calculated from single-cell dataset sampled at E25 (S9).(f) TPMs from the bulk RNAseq for MLANA and GPNMB, both used to perform whole-mount in situ hybridization experiments.The elements of the box plots are: centre line, median; box limits, upper and lower quartiles; whiskers, 1.5x interquartile range; points, measurements).Source data are provided as a Source Data file.
results of the PMEL target regions for functional analyses.Electropherogram of the target region sequenced on DNA extracted from (a) control corn snake fibroblasts that were transfected with a blank (no gRNA), (b) corn snake fibroblasts transfected with gRNA 204f, (c) corn snake fibroblasts transfected with gRNA 302r, (d) corn snake PMEL mutant with a wild type phenotype (204f), (e) corn snake PMEL mutant with stripes (204f), (f) corn snake PMEL mutant with both stripes and blotches (204f), (g) biopsy 1 taken from (f) from a region with stripes, (h) biopsy 2 taken from (f) from a region with blotches, (i) corn snake PMEL mutant with a wild type phenotype (302r), (j) corn snake PMEL mutant with stripes (302r), (k) corn snake PMEL mutant with a wild type phenotype (302r) and heterozygous, (l) cloned allele 1 from (k), (m) cloned allele 2 from (k), (n) corn snake PMEL mutant with a modified phenotype (302r), (o) biopsy 1 taken from (n) from an anterior region, (p) biopsy 2 taken from (n) from a posterior region.WT: wild type, S: stripes, B: blotches, red rectangles: PAM site, red arrows: insertion/deletion sites.(collagen type I alpha 1 chain) and COL1A2 (collagen type I alpha 2 chain).(d) Count of genes that intersect between the top 100 genes expressed in the two fibroblast clusters and the gene list specific to different mouse fibroblast types, as defined by Jacob et al., 2023.(e) UMAP embedding showing the expression of fibroblast markers in the two fibroblast clusters.Supplementary Figure 11.Xanthophores and iridophores markers.(a) In situ hybridization of a PAX7 probe, a classical marker of xanthophores in zebrafish, on a S9 wildtype embryo.With this method, we could detect very faint expression which seems to be in the background region, and not within the blotches.The regions that are likely to correspond to the blotches are delineated by a white dashed line.This result should be taken with caution.Species-specific xanthophores markers are needed for confirmation.(b) Expression levels of classical xanthophore (CSF1R, GCH1) and iridophore markers (TFEC, CSF1) are very low at this stage of development in corn snakes.

Supplementary Table 1. Genomic DNA libraries used for mapping the Terrazzo causative variant.
In parentheses, we provide the total number of individuals in each library (column 'Source') and the average coverage for a 1.7 Gb genome (column 'Filtered reads').PE: Pairedend reads.

Table 2 . List of the 35 protein coding genes in the Terrazzo interval on Super-scaffold 85.
We provide information on the presence of co-segregating variants in the coding region of each gene ('SNPs/MNPs/indels' column) and of any non-synonymous aminoacid modifications in the protein coding part of genes detected in the Terrazzo libraries when compared to the wild-type ('aa modifications' column).We include information on APOH, a differentially expressed gene in Terrazzo animals and located near, but not within, the reduced genomic interval.