Stripes and loss of color in ball pythons (Python regius) are associated with variants affecting endothelin signaling

Abstract Color patterns in nonavian reptiles are beautifully diverse, but little is known about the genetics and development of these patterns. Here, we investigated color patterning in pet ball pythons (Python regius), which have been bred to show color phenotypes that differ dramatically from the wildtype form. We report that several color phenotypes in pet animals are associated with putative loss-of-function variants in the gene encoding endothelin receptor EDNRB2: (1) frameshift variants in EDNRB2 are associated with conversion of the normal mottled color pattern to skin that is almost fully white, (2) missense variants affecting conserved sites of the EDNRB2 protein are associated with dorsal, longitudinal stripes, and (3) substitutions at EDNRB2 splice donors are associated with subtle changes in patterning compared to wildtype. We propose that these phenotypes are caused by loss of specialized color cells (chromatophores), with loss ranging from severe (fully white) to moderate (dorsal striping) to mild (subtle changes in patterning). Our study is the first to describe variants affecting endothelin signaling in a nonavian reptile and suggests that reductions in endothelin signaling in ball pythons can produce a variety of color phenotypes, depending on the degree of color cell loss.


Introduction
Skin colors in vertebrates originate from specialized color cells.These cells are known as chromatophores, and they come in three main types (Kelsh et al. 2009;Mills and Patterson 2009;Cuthill et al. 2017).Melanophores synthesize the pigment melanin, whose color ranges from black to reddish brown.Xanthophores store red-to-yellow pigments, including pteridines synthesized by the cell and carotenoids obtained from the diet.Iridophores contain platelets of crystallized purines, which can appear blue, white, or iridescent, depending on their structure.Skin color in nonavian reptiles and lower vertebrates is determined by layering of different chromatophore types (Patterson and Parichy 2019;Kuriyama et al. 2020).Green colors, for example, are produced by yellow xanthophores above blue iridophores (e.g.Saenko et al. 2013).Golden-brown colors, as a second example, are produced by a combination of brown melanophores, yellow xanthophores, and white iridophores (e.g.Brejcha et al. 2019).Mammals possess only one type of chromatophore (melanophores), and their skin color is limited to shades of brown and black (Mort et al. 2015).
Chromatophores are derived from the neural crest and migrate during development to reach their final positions in the skin (Mort et al. 2015;Parichy and Spiewak 2015;Haupaix and Manceau 2020).Once in the skin, chromatophores organize into patterned arrangements to produce the overall color pattern of the animal.This process has been best characterized in zebrafish, where patterning is driven by short-and long-range interactions among chromatophores themselves (Singh and Nuesslein-Volhard 2015;Patterson and Parichy 2019).Similar mechanisms likely function in mammals to produce periodic markings such as stripes and spots (Kaelin et al. 2012(Kaelin et al. , 2021;;Mallarino et al. 2016).Pattern formation is less well understood in nonavian reptiles, although the idea that similar mechanisms may function in this group is supported by theoretical studies showing that Turing-like interactions can produce many of the basic color patterns observed in this group (Murray and Myerscough 1991;Kondo and Shirota 2009;Kuriyama et al. 2020;Ascarrunz and Sanchez-Villagra 2022).
An emerging model for understanding coloration and pattern formation in nonavian reptiles is the ball python (Python regius) (Garcia-Elfring et al. 2023;Brown et al. 2022).Ball pythons are native to west Africa, but have become common as pets in the United States.Wild ball pythons exhibit a mottled color pattern, consisting of irregular black and golden-brown patches of skin (Fig. 1).Pet ball pythons, by contrast, have been bred to include animals with dramatically different color patterns (McCurley 2005;Broghammer 2019).These variants, known as "color morphs," include animals with solid body colors, longitudinal stripes, or major changes in pattern patchiness.Many of these phenotypes are reminiscent of patterning changes that have occurred within the snake lineage across evolutionary time, including patterns thought to be adaptive for certain behavioral ecologies (Brodie 1992;Allen et al. 2013;Olsson et al. 2013;Pizzigalli et al. 2020).
Color morphs in ball pythons may therefore provide insight into the genetic and developmental mechanisms controlling patterning changes across species.
Many color morphs in ball pythons show simple dominant or recessive patterns of inheritance.These inheritance patterns suggest that these color morphs are caused by variants in single genes (e.g.Garcia-Elfring et al. 2023;Brown et al. 2022).Some morphs belong to allelic series, with some morphs in the series showing more dramatic phenotypes than others (McCurley 2005;Broghammer 2019).Such morphs are thought to represent different molecular variants of the same gene, akin to mutants in allelic series arising from classical mutagenesis screens in model organisms.
One example of an allelic series in ball pythons is the Yellowbelly series of color morphs (Fig. 1).This series is described by breeders as having five alleles (yellowbelly, spark, specter, gravel, and asphalt).Morphs in the Yellowbelly series exhibit changes in patterning and reduced coloration compared to wildtype.These phenotypes range from severe to mild (Fig. 1a and b).Homozygotes of the strongest allele (yellowbelly) have skin that is almost fully white (Fig. 1b).The intermediate alleles (spark and specter) produce homozygotes with a pair of dorsal, longitudinal stripes (Fig. 1b).These animals also have moderately reduced coloration.The weakest alleles (gravel and asphalt) produce homozygotes with subtle alterations in patterning and mildly reduced coloration (Fig. 1b).Compound heterozygotes show phenotypes intermediate between the parents, often including dorsal, longitudinal stripes (Fig. 1d).These phenotypes are largely recessive, and heterozygotes of any of these alleles in combination with the wildtype allele show color patterns that are largely similar to wildtype (Fig. 1f).Some heterozygotes have slightly lighter underbellies than wildtype (Fig. 1e and f), but this phenotype is variable and difficult to recognize, especially for alleles other than the yellowbelly allele.
All-white skin in ball pythons is not the outcome of a defect in pigment synthesis.Such defects remove a single pigment (e.g.melanin) without affecting other pigments (e.g. Brown et al. 2022).Instead, all-white skin suggests a loss of multiple types of chromatophores.The absence of brown-to-black coloration suggests an absence of melanophores, and the absence of yellow coloration suggests an absence of xanthophores.Iridophores in animals with all-white skin are harder to asses, given that iridophores in ball pythons are likely white-reflecting and thus cannot be detected by visual appearance alone.In all-white morphs of other reptiles and amphibians, for example, iridophores are present in some cases (e.g.Ullate-Agote and Tzika 2021) and absent in others (e.g.Woodcock et al. 2017).Loss of chromatophores has been described in other vertebrates, where the underlying defects include lack of chromatophore specification, insufficient proliferation of chromatophore precursors, and failure of these precursors to migrate from the neural crest to their normal positions in the skin (reviewed in Mort et al. 2015;Patterson and Parichy 2019).
A prominent phenotype in the Yellowbelly series is conversion of the wildtype color pattern to dorsal, longitudinal stripes (Fig. 1b).We hypothesized that this phenotype, like all-white skin, was the outcome of chromatophore loss, albeit to a lesser degree.This hypothesis arises from a resemblance between dorsal striping in ball pythons and a mutant phenotype in zebrafish caused by fewer chromatophores.The adult color pattern in zebrafish is initiated by iridophore precursors that migrate via the horizontal midline to populate the skin (Fig. 2a) (Singh and Nuesslein-Volhard 2015;Patterson and Parichy 2019).In mutant zebrafish with fewer iridophores, iridophores fail to populate the ).Xanthophore precursors persist in the epidermis from larvalhood.Melanophore and iridophore precursors migrate from the neural crest to the epidermis along nerve tracks.The primary route for iridophores is through the horizontal midline.The adult pattern consists of dark stripes (melanophores and sparse iridophores) and light interstripes (xanthophores and dense iridophores).b) Schematic of mutant zebrafish with a primary loss of iridophores and a secondary loss of melanophores (Parichy et al. 2000;Lopes et al. 2008;Frohnhöfer et al. 2013;Patterson and Parichy 2013;Mo et al. 2017;Spiewak et al. 2018).Iridophores fail to populate the skin and remain clustered near the horizontal midline.Melanophores require iridophores for survival and persist only near the iridophores.c) Hypothesized model for chromatophore migration in ball pythons.Chromatophores are presumed to migrate in a dorsal-to-ventral path, as they do in other higher vertebrates (Mort et al. 2015).d-f)Phenotypes in the Yellowbelly series are hypothesized to be caused by a loss of chromatophores.Loss is hypothesized to range from mild (altered patterning) to moderate (striped) to severe (white).This model depicts a loss of all three types of chromatophores, but the loss of iridophores is uncertain.The range of developmental processes plausibly affected includes chromatophore specification, proliferation, migration, or a combination thereof.U. M. Dao et al. | 3 skin and instead remain clustered at the horizontal midline (Parichy et al. 2000;Lopes et al. 2008;Patterson and Parichy 2013;Frohnhöfer et al. 2013;Mo et al. 2017;Spiewak et al. 2018).Melanophores in zebrafish require iridophores for survival, and hence melanophores in these mutants are less abundant and persist only near iridophores (Frohnhöfer et al. 2013;Patterson and Parichy 2013;Patterson et al. 2014).The resulting phenotype is a pair of midline stripes on each side of the fish's body (Fig. 2b).We suspected that these midline stripes might be analogous to dorsal stripes in the Yellowbelly series, given that chromatophores in ball pythons likely migrate from the neural crest in a dorsal-to-ventral path, as they do in other higher vertebrates (Mort et al. 2015).Incomplete loss of chromatophores in ball pythons might therefore cause chromatophores to remain clustered near the dorsal ridge.This hypothesis provides a simple explanation for the range of phenotypes observed in the Yellowbelly series: all-white skin may be outcome of chromatophore loss that is severe (Fig. 2f); dorsal striping may be the outcome of chromatophore loss that is moderate (Fig. 2e); and subtle alternations in patterning may be the outcome of chromatophore loss that is mild (Fig. 2d).The cellular basis of this loss could plausibly include defects in chromatophore specification, proliferation, or migration, or a combination thereof.
Our hypothesis about chromatophore loss in the Yellowbelly series suggested candidate genes that might be causative for these phenotypes.At the top of this list were genes encoding components of endothelin signaling.Endothelin signaling comprises a vertebrate-specific signaling pathway responsible for regulating the development of neural crest cells (Saldana-Caboverde and Kos 2010; Bondurand et al. 2018).Components of endothelin signaling include endothelin ligands (endothelins), G-proteincoupled receptors (endothelin receptors), and enzymes responsible for processing the ligands (endothelin-converting enzymes).Most vertebrate genomes encode three endothelin ligands (EDN1, EDN2, and EDN3) and three endothelin receptors (EDNRA, EDNRB1, and EDNRB2) (Braasch et al. 2009).Exceptions include zebrafish and eutherian mammals, which have lost EDNRB2 (Braasch et al. 2009).EDNRA is selective for ligands EDN1 and EDN2 and helps pattern the skeleton of the face and head (Krystek et al. 1994;Lee et al. 1994;Lecoin et al. 1998;Liu et al. 2019; but see Menzi et al. 2016;Clouthier et al. 2010).EDNRB1 and EDNRB2, by contrast, regulate chromatophore precursors, among other cell types.EDNRB1 is required for proper proliferation and migration of chromatophore precursors in zebrafish and mammals, and EDNRB2 plays a similar role in birds (e.g.Baynash et al. 1994;Parichy et al. 2000;Pla et al. 2005;Harris et al. 2008;reviewed in Braasch and Schartl 2014).EDNRB1 and EDNRB2 are activated in cell culture by all three endothelin ligands (Krystek et al. 1994;Lecoin et al. 1998;Liu et al. 2019), but their primary ligand in vivo is thought to be EDN3.Loss of EDN3 phenocopies loss of EDNRB1 in zebrafish and mammals (Baynash et al. 1994;Edery et al. 1996;Hofstra et al. 1996;Bondurand et al. 2006;Pingault et al. 2010;Spiewak et al. 2018), and reduced EDN3 expression causes similar phenotypes in amphibians (Kawasaki-Nishihara et al. 2011;Woodcock et al. 2017).EDN3 is a potent mitogen of chromatophore precursors across vertebrates (Lahav et al. 1996;Reid et al. 1996;Dupin et al. 2000), and increased EDN3 expression is associated with darker skin and hair in birds and mammals (Dorshorst et al. 2011;Shinomiya et al. 2012;Kaelin et al. 2021).EDNRB1, EDNRB2, and EDN3 therefore play conserved roles in chromatophore development in vertebrates and are candidate genes for the causative gene of the Yellowbelly series.
The goal of the current study was to identify the genetic cause of color phenotypes in the Yellowbelly series.We predicted that these phenotypes were caused by loss-of-function variants in the endothelin ligand gene EDN3 or the endothelin receptor genes EDNRB1 or EDNRB2.Our approach was to recruit samples of pet ball pythons (shed skins) from across the United States and search for putative loss-of-function variants in these genes.

Recruitment of ball python sheds
Ball python sheds were recruited from pet owners and breeders by placing announcements in Twitter, Reddit, Instagram, and Facebook, and by contacting sellers having active listings on Morph Market (www.morphmarket.com).Contributors were instructed to allow sheds to dry (if wet) and to ship sheds via standard first-class mail.Contributors sending multiple sheds were instructed to package sheds individually in plastic bags during shipping.Contributors were not provided monetary compensation for sheds, although some contributors were given prepaid shipping envelopes to cover shipping costs.Upon receipt, sheds were stored at −20°C to kill any insect larvae infesting the sheds.
We attempted to maximize genetic diversity within each category of morph by excluding animals described as siblings, parents, or offspring of another animal in our sample.However, we cannot exclude the possibility that some animals in our sample may have been close relatives, given that full pedigree information was not available for most animals.The idea that animals in our sample were derived from multiple lineages is supported by our discovery of multiple molecular variants corresponding to the allele known to breeders as yellowbelly.Owner codes and US states of residence for each animal are provided in Supplementary Table 2.
The total set of animals used in our study were 70 animals described as yellowbelly homozygotes, 130 animals described as yellowbelly heterozygotes, 22 animals described as spark/yellowbelly, 3 animals described as specter homozygotes, 27 animals described as specter/yellowbelly, 7 animals described as gravel homozygotes, 22 animals described as gravel/yellowbelly, 4 animals described as asphalt homozygotes, 23 animals described as asphalt/yellowbelly, and 48 animals described as having wildtype coloration or belonging to morphs other than those in the Yellowbelly series (Non-Yellowbelly animals).

Performing experiments in an undergraduate laboratory course
The majority of experiments and analyses described in this study were performed by undergraduate students as part of a laboratory course at Eastern Michigan University (BIO306W).This practice required that our experimental design rely on simple techniques, namely PCR, restriction digests, and Sanger sequencing.To avoid student errors in these techniques, we implemented the following precautions.First, students never handled more than one shed skin at the same time.Second, students performed negative and positive control reactions for all PCR amplifications and restriction digests.Data from students having incorrect controls were excluded from analysis.Third, all sequence analyses were performed independently by three or more students.When the results of all students did not all agree, sequences were re-analyzed by the instructor (H.S.S.).Fourth, when we encountered an animal whose molecular genotype did not match the breeder description (e.g. the single animal labeled spark/yellowbelly whose genotype matched specter/yellowbelly and vice versa-see Results and Discussion), DNA from the animal was re-extracted by the instructor (H.S.S.), and the molecular genotypes were determined a second time.This precaution was taken to guard against mix-ups by students, but it proved unnecessary because the second round of genotypes always matched the first.

DNA extraction
Sheds containing visible dirt or debris were rinsed in tap water until clean.Sheds were air-dried and lysed overnight at ∼60°C in ∼1 ml lysis buffer (100 mM Tris-HCl pH 8.0, 100 mM EDTA, 2% sodium dodecyl sulfate, 3 mM CaCl 2 , 2 mg/ml Proteinase K) per palm-sized piece of shed.Lysis reactions were performed in standard 1.5 ml tubes.In some cases, pieces of shed were snipped or torn into multiple pieces to enable them to fit into tubes more easily.Dorsal and ventral scales were used interchangeably, after preliminary studies indicated that DNA could be reliably extracted from both types of scales.
Lysate was separated from visible fragments of undigested shed and further cleared by centrifugation at 13,000×g for 2 min.To precipitate protein, ammonium acetate was added to supernatant to a final concentration of 1.875-2.5M. Samples were incubated on ice for 5 min and centrifuged at 13,000×g for 3 min at 4°C.Supernatant was mixed with an equal volume of magnetic-bead mixture [10 mM Tris-HCl pH 8.0, 1 mM EDTA, 1.6 M NaCl, 0.2% Tween-20, 11% polyethylene glycol, 0.04% washed SpeedBeads (Sigma #GE45152105050250)] and shaken for 5 min.Beads were separated from supernatant using a magnet, washed twice in 0.2 ml 70% ethanol for 2 min, and air-dried for ∼1 min.DNA was eluted from beads in TE buffer (10 mM Tris-HCl pH 8.0, 1 mM EDTA) at 65°C for >5 min.

Primer design and PCR
Primers were designed against the genome of Burmese python, the closest relative of ball python for which genome sequence was available (Castoe et al. 2013).Annotations for Burmese python were XM_025169570 (mRNA) and XP_007432959 (protein) for EDN3, XM_007424069 (mRNA) and XP_007424131 (protein) for EDNRB1, and XM_007436628 (mRNA) and XP_007436690 (protein) for EDNRB2.These annotations were manually examined for errors, by compared gene structure across vertebrates.For both genes, annotations were determined to be correct, based on the observation that gene structure matched gene structure in other vertebrates, including in species with higher quality annotations, such as mouse, chicken, and zebrafish.
Primers were designed using Primer3 (Untergasser et al. 2012), using default parameters and a target annealing temperature of 60°C.Amplification was performed using an annealing temperature of 57°C, to allow for occasional divergence between ball python and Burmese python genomic sequences.Coding regions of EDNRB1, EDNRB2, and EDN3 were amplified and sequenced using primers given in Supplementary Table 1.

The yellowbelly allele is associated with frameshift variants in EDNRB2
We hypothesized that morphs in the Yellowbelly series were caused by loss-of-function variants in EDN3, EDNRB1, or EDNRB2.To test this hypothesis, we amplified and sequenced the coding regions and adjacent splice sites of each gene in one ball python described as a yellowbelly homozygote by its owner.Sequences from this animal were compared to sequences from an animal having normal coloration (henceforth "wildtype").We observed that the yellowbelly homozygote was homozygous for a 1-bp deletion in the fourth coding region of EDNRB2 (OP589186: c.1646del, Fig. 3b).No coding or splice-site variants were observed in EDN3 or EDNRB1.The 1-bp deletion in EDNRB2 introduces a frameshift into the transcript, which removes or alters three of the seven transmembrane helices of the EDNRB2 protein (Fig. 3e).This truncation likely abolishes protein function, given that G-protein-coupled receptors are generally unable to function without an intact transmembrane structure (Rosenbaum et al. 2009).We conclude that the 1-bp deletion is likely a strong loss-of-function allele of EDNRB2.
To test whether the 1-bp deletion in EDNRB2 was associated with phenotypes attributed to the yellowbelly allele, we genotyped 117 additional animals for this deletion.These animals consisted of 69 animals described as yellowbelly homozygotes and 48 animals described as having normal coloration or belonging to morphs other than those in the Yellowbelly series (henceforth "Non-Yellowbelly" animals).These animals were collected from a total of 35 and 33 pet owners, respectively, and thus represent a broad sampling of the pet population.We observed that all animals described as yellowbelly homozygotes were either homozygous for the 1-bp deletion (60 animals) or heterozygous for it (nine animals) (Fig. 3c).By contrast, none of the Non-Yellowbelly animals carried the deletion (Fig. 3c).These findings demonstrate a strong but imperfect association between the 1-bp deletion and morph type.
We hypothesized that the nine animals described as yellowbelly homozygotes but heterozygous for the 1-bp deletion might carry an alternate loss-of-function variant in EDNRB2 on the homologous chromosome.To test this hypothesis, we sequenced the coding regions and adjacent splice sites of EDNRB2 in one of these animals.We observed that this animal was heterozygous for a 17-bp deletion in the fourth coding region of EDNRB2 (OP589186:c.1747_1763del,Fig. 3b).This deletion was located on the homologous chromosome compared to the 1-bp deletion (Supplementary Fig. 1).The 17-bp deletion introduces a frameshift into the transcript, which removes or alters two of the seven transmembrane helices of the EDNRB2 protein (Fig. 3f).The 17-bp deletion thus represents a second putative loss-of-function allele of EDNRB2.
To assess the contribution of the 17-bp deletion to phenotypes attributed to the yellowbelly allele, we genotyped all animals in our panel for the 17-bp deletion.We found that the 17-bp deletion was present exclusively in the nine yellowbelly homozygotes heterozygous for the 1-bp deletion (Fig. 3c).In each case, the 17-bp deletion was carried on the homologous chromosome compared to the 1-bp deletion (Supplementary Table 2).These results show that every animal described as a yellowbelly  showing association between the yellowbelly allele and the 1-bp and 17-bp deletions.This table includes the original yellowbelly homozygote in which the 1-bp deletion was identified.Numbers of owners and US states indicate the total numbers of owners and US states from which the animals were collected.*, presumed mislabeling by owner, including one animal heterozygous for L152F (see Discussion).d) Schematic of wildtype EDNRB2 protein, created using Protter (Omasits et al. 2014).e, f) Schematics of truncated EDNRB2 proteins resulting from the 1-bp and 17-bp deletions.Schematics are intended to show premature terminations and do not imply that the truncated proteins fold properly and are inserted into the cell membrane.
homozygote was either homozygous for the 1-bp deletion or compound heterozygous for the 1-bp deletion and the 17-bp deletion (Fig. 3d).The absence of animals homozygous for the 17-bp deletion may simply reflect the low frequency of this allele in our sample.These data suggest that the allele known to breeders as yellowbelly represents the collection of two molecular alleles of EDNRB2: a 1-bp deletion and a 17-bp deletion.We propose that these deletions produce the same phenotype and were therefore not recognized by breeders as distinct.

The spark and specter alleles are associated with EDNRB2 missense variants
Animals described as carrying the spark or specter alleles of the Yellowbelly series show a more moderate loss of coloration than animals described as yellowbelly homozygotes (Fig. 1b).We therefore predicted that the spark and specter alleles might represent variants in EDNRB2 that reduce gene function to a lesser degree than do the 1-bp and 17-bp deletions.To test this hypothesis, we sequenced the coding regions and adjacent splice sites of EDNRB2 in two animals: one animal described as a spark/yellowbelly compound heterozygote and one animal described as a specter homozygote.(Our sample did not include any animals described as spark homozygotes.)We observed that the spark/ yellowbelly animal was heterozygous for a missense variant in the first coding region of EDNRB2 (OP589186:c.G481C).This variant leads to a leucine-to-phenylalanine exchange at position 152 of the protein, termed hereafter L152F (Fig. 4c).This animal was also heterozygous for the 1-bp deletion associated with the  showing associations between the spark and specter alleles and the missense variants.This table includes the original spark/yellowbelly compound heterozygote and the original specter homozygote in which the missense variants were identified.Numbers of owners and US states indicate the total numbers of owners and US states from which the animals were collected.*, presumed mislabeling by owner (see Discussion).e, f) Schematics of EDNRB2 proteins containing missense variants, created using Protter (Omasits et al. 2014).g) Alignment of EDNRB2 and EDNRA protein sequences surrounding missense variants.EDNRA is included to demonstrate conservation outside EDNRB2.Full EDNRA sequences were not available for ball python or painted turtle.
yellowbelly allele, as expected.We observed that the specter homozygote was homozygous for a missense variant in the fifth coding region of EDNRB2 (OP589186:c.G2601C).This variant leads to an arginine-to-proline exchange at position 315 of the protein, termed hereafter R315P (Fig. 4c).No other coding or splice-site variants were observed in either animal.Both missense variants affected conserved sites of the EDNRB2 protein, suggesting that they likely disrupt protein function.Position 152, which is located in a transmembrane helix (Fig. 4e), is conserved as a leucine or isoleucine across species (Fig. 4g).Position 315, which is located in an intracellular loop (Fig. 4f), is conserved as an arginine across species (Fig. 4g).We conclude that L152F and R315P are candidates for the spark and specter alleles, respectively.
To test whether L152F and R315P were associated with phenotypes attributed to the spark and specter alleles, we genotyped these variants in the following panel of animals: 21 additional animals described as spark/yellowbelly compound heterozygotes; two additional animals described as specter homozygotes; 27 animals described as specter/yellowbelly compound heterozygotes; and 48 Non-Yellowbelly animals.Animals were also genotyped for the 1-bp and 17-bp deletions associated with the yellowbelly allele.We observed that the missense variants were nearly perfectly associated with the spark and specter alleles: (1) both animals described as specter homozygotes were homozygous for R315P, (2) all but one animal described as spark/yellowbelly were heterozygous for L152F and for one of the deletions, and (3) all but one animal described as specter/yellowbelly were heterozygous for R315P and for one of the deletions (Fig. 4d).The exceptions were that one animal described as spark/yellowbelly carried R315P, rather than L152F, and one animal described as specter/yellowbelly carried L152F, rather than R315P (Fig. 4d).These exceptions may represent mislabeling by owners, given the phenotypic similarity between the spark and specter alleles (see Discussion).These data largely support the hypothesis that L152F represents the spark allele and R315P represents the specter allele.

The gravel and asphalt alleles are associated with variants in EDNRB2 splice donors
Animals described as carrying the gravel or asphalt alleles of the Yellowbelly series show the mildest loss of coloration of all morphs in the series (Fig. 1b).We therefore predicted that the gravel and asphalt alleles might represent variants in EDNRB2 with mild effects on gene function.To test this hypothesis, we sequenced the coding regions and adjacent splice sites of EDNRB2 in one animal described as a gravel homozygote and one animal described as an asphalt homozygote.We observed that both animals were homozygous for synonymous G-to-A nucleotide substitutions in the −1 positions of splice donors.For the gravel homozygote, the substitution was located in the splice donor for the first intron (OP589186:c.G499A, Fig. 5c).For the asphalt homozygote, the substitution was located in the splice donor for the sixth intron (OP589186:c.G3118A, Fig. 5c).No coding variants and no other splice-site variants were observed in either animal.The −1 position of splice donors contributes to splicing efficiency across eukaryotes (Cartegni et al. 2002), and this position is a guanine in ∼80% of splice donors in vertebrate genes (Ma et al. 2015).In EDNRB2, the −1 positions of the first and sixth splice donors are conserved as guanines across species (Fig. 5c), suggesting that guanines at these sites contribute to gene function.We conclude that the splice-donor substitutions identified herein are candidates for the gravel and asphalt alleles.
To test whether the G-to-A substitutions in splice donors were associated with phenotypes attributed to the gravel and asphalt alleles, we genotyped these variants in the following panel of animals: six additional animals described as gravel homozygotes; three additional animals described as asphalt homozygotes; 22 animals described as gravel/yellowbelly compound heterozygotes; 23 animals described as asphalt/yellowbelly compound heterozygotes; and 48 Non-Yellowbelly animals.These animals were also genotyped for the 1-bp and 17-bp deletions associated with the yellowbelly allele.We observed that the G-to-A substitutions in splice donors were nearly perfectly associated with the gravel and asphalt alleles: (1) all animals described as gravel or asphalt homozygotes were homozygous for the substitutions in the first or sixth splice donors, respectively, (2) all animals described as gravel/yellowbelly were heterozygous for the substitution in the first splice donor and for one of the deletions, and (3) all but two animals described as asphalt/yellowbelly were heterozygous for the substitution in the sixth splice donor and for one of the deletions (Fig. 5d).The exceptions were two animals described as asphalt/yellowbelly, which carried the substitution in the first splice donor, rather than in the sixth splice donor (Fig. 5d).These exceptions may represent mislabeling by owners, given the phenotypic similarity between the gravel and asphalt alleles (see Discussion).These data largely support the hypothesis that the G-to-A substitution in the first splice donor represents the gravel allele, and the G-to-A substitution in the sixth splice donor represents the asphalt allele.

No candidate variants in EDN3 or EDNRB1
Our initial focus on EDNRB2, rather than EDN3 and EDNRB1, relied on sequences from a single yellowbelly homozygote.To further exclude EDN3 and EDNRB1 as causative for phenotypes in the Yellowbelly series, we amplified and sequenced the coding regions and adjacent splice sites of EDN3 and EDNRB1 in the following animals: One specter homozygote (R315P/R315P), one gravel homozygote (A/A at the first splice donor), one asphalt homozygote (A/A at the sixth splice donor), and one spark/yellowbelly compound heterozygote (L152F/1-bp deletion).We observed no coding variants and no splice-site variants in any of these animals.These results further support the conclusion that phenotypes in the Yellowbelly series are caused by variants in EDNRB2, not EDN3 or EDNRB1.

Reliability of heterozygote descriptions
Ball pythons described as heterozygous for alleles in the Yellowbelly series sell for higher prices than wildtype animals.This price difference provides an incentive for sellers to label animals as heterozygotes.Yet heterozygotes for alleles in the Yellowbelly series are difficult to distinguish from wildtype.Heterozygotes have subtly lighter coloration and altered patterning in the underbelly (Fig. 1e and f), but this phenotype is variable and overlaps with the range of phenotypes observed among wildtype animals.Thus, animals described as heterozygous for alleles in the Yellowbelly series are viewed with a degree of skepticism.Many believe that some of these animals are mislabeled, perhaps for monetary gain by sellers.
Our discovery of association between alleles in the Yellowbelly series and variants in EDNRB2 allowed us to test the reliability of heterozygote descriptions.To do so, we genotyped 130 animals described as yellowbelly heterozygotes for the 1-bp and 17-bp deletions in EDNRB2.These animals were collected from a total of 68 pet owners across the United States and therefore represent a broad sampling of the pet population.We focused on the yellowbelly allele, rather than other alleles in the Yellowbelly series, because the yellowbelly allele is described as producing the strongest phenotype in heterozygotes.We observed that 98 of the 130 animals in our sample were heterozygous for the 1-bp or 17-bp deletions in EDNRB2 (Fig. 3c).The animals lacking these deletions were further tested to determine whether they carried any of the other variants associated with alleles in the Yellowbelly series-L152F (spark), R315P (specter), and the G-to-A substitutions in splice donors (gravel and asphalt).One animal was found to be heterozygous for L152F (Supplementary Table 2).The other animals did not carry any of these variants (Supplementary Table 2).Thus, of the 130 animals described as yellowbelly heterozygotes, 98 were correctly described, one was a spark heterozygote, and the remainder carried only putatively wildtype copies of EDNRB2.These findings confirm the skepticism surrounding heterozygote descriptions for alleles in the Yellowbelly series.They suggest that approximately one in four animals described as a yellowbelly heterozygote is mislabeled.

Discussion
The goal of the current study was to identify the genetic cause of phenotypes in the Yellowbelly series of ball pythons.These phenotypes were characterized by a loss of coloration and changes in patterning, ranging from severe (all-white skin) to moderate (dorsal striping) to mild (subtle changes in patterning compared to wildtype) (Fig. 1).Through careful choice of candidate genes and a bit of luck, we discovered that these phenotypes were associated with six putative loss-of-function variants in the gene encoding endothelin receptor EDNRB2 (Fig. 6).All-white skin was associated with frameshift variants (Fig. 3).Dorsal striping was associated with missense variants (Fig. 4).Subtle changes in patterning were associated with splice-donor variants (Fig. 5).To our knowledge, these variants represent the first description of genetic variants affecting endothelin signaling in a nonavian reptile.

Causality, limitations, and the candidate-gene approach
We propose that the EDNRB2 variants described herein are causal for phenotypes in the Yellowbelly series.This idea is supported by a correlation between the severity of phenotypes in the Yellowbelly series and the severity of genetic lesions in EDNRB2  showing associations between the gravel and asphalt alleles and the nucleotide substitutions in splice donors.This table includes the original gravel and asphalt homozygotes in which the nucleotide substitutions were identified.Numbers of owners and US states indicate the total numbers of owners and US states from which the animals were collected.*, presumed mislabeling by owner (see Discussion).e) Alignment of nucleotide sequences for the first and sixth splice donors of EDNRB2.
(Fig. 6).The most severe phenotype in the series-all-white skinwas associated with frameshift variants in EDNRB2, which plausibly represent null alleles (Fig. 3).White skin is characteristic of an absence of melanophores and xanthophores (e.g.Woodcock et al. 2017; Ullate-Agote and Tzika 2021), and the association between all-white skin and EDNRB2 in ball pythons is consistent with EDNRB2 controlling chromatophore development throughout vertebrates (reviewed in Braasch and Schartl 2014).We propose that all-white skin represents the null phenotype for EDNRB2 in ball pythons.
The next most severe phenotype in the Yellowbelly series was dorsal striping.This phenotype was associated with missense variants at conserved sites of the EDNRB2 protein (Fig. 4).We propose that these variants represent hypomorphic alleles.These missense variants might reduce EDNRB2 protein function via effects on folding or interactions with cofactors, among other possibilities.Dorsal striping in ball pythons resembles mutant phenotypes in zebrafish caused by iridophores and melanophores being reduced in number and remaining clustered near the horizontal midline (Fig. 2a and b) (Parichy et al. 2000;Lopes et al. 2008;Frohnhöfer et al. 2013;Patterson and Parichy 2013;Mo et al. 2017;Spiewak et al. 2018).We propose that stripe formation in ball pythons is analogous to this phenotype in zebrafish, and that stripes in ball pythons are caused by chromatophores remaining clustered near the dorsal ridge (Fig. 2e).This phenotype might arise from reduced cell proliferation, reduced cell migration, or both, given that endothelin signaling controls one or both processes in other species (Lahav et al. 1996;Reid et al. 1996 The mildest phenotype in the Yellowbelly series was altered patterning.This phenotype was associated with G-to-A nucleotide substitutions in the −1 positions of EDNRB2 splice donors (Fig. 5).The −1 position of splice donors contributes to splicing efficiency (Cartegni et al. 2002), and variants in the −1 positions of splice donors of human genes have been known to disrupt gene function (Anna and Monika 2018).We propose that the splice-donor substitutions in EDNRB2 reduce the fidelity of splicing and thereby reduce EDNRB2 protein levels.The outcome may be a mild reduction in chromatophore number or altered cell migration or both, leading to changes in patterning (Fig. 2d).The idea that mild changes cell number or migration could produce noticeable changes in patterning is consistent with the idea that patterns arise via a complex network of cell-to-cell interactions (Singh and Nuesslein-Volhard 2015;Patterson and Parichy 2019), which are sensitive to cell numbers (Spiewak et al. 2018) and other parameters (Kaelin et al. 2012(Kaelin et al. , 2021;;Watanabe and Kondo 2012;McGuirl et al. 2020).
The scope of our study was limited to three genes: EDN3, EDNRB1, and EDNRB2.Thus, we cannot formally exclude the possibility that variants in other regions of the genome contribute to phenotypes in the Yellowbelly series.However, we view this possibility as unlikely, given that six independent haplotypes of EDNRB2 were associated with phenotypes in the Yellowbelly series, and that the severity of these phenotypes correlated with the severity of genetic lesions in EDNRB2.We propose that further evidence supporting a link between EDNRB2 and the Yellowbelly series could be obtained by introducing the variants described herein into a model system (e.g.zebrafish) to test for effects on protein function or gene expression.An additional limitation of our study is that samples were limited to shed skins.Shed skins did not allow us to compare chromatophore numbers across morphs, nor did they allow us to characterize chromatophore behavior during development.Characterization of embryos from wildtype ball pythons, as well as embryos in the Yellowbelly series, is needed to fully characterize the role of EDNRB2 in ball pythons.

Pros and cons of community-sourced sampling
Our study relied on pet samples recruited from the community, rather than samples from a single pedigree.An advantage of this approach was increased genetic diversity in our sample.An example of this diversity was the discovery that the allele known to breeders as yellowbelly represents two distinct molecular alleles of EDNRB2-a 1-bp deletion and a 17-bp deletion (Fig. 3).The 17-bp deletion was approximately sixteen-fold less common in our sample than the 1-bp deletion and might have been plausibly missed by our study, had we focused on a single pedigree.The phenomenon of multiple molecular alleles corresponding to a single breeder-defined allele is not unique to the Yellowbelly series and has been reported for two other color morphs in ball pythons-the Albino color morph and the Ultramel color morph (Brown et al. 2022).This phenomenon highlights the utility of community-sourced sampling in ball pythons and shows that the pet population is more diverse than previously recognized.
A disadvantage of using pet samples sourced from the community is that genotype descriptions provided by owners come with a degree of uncertainty.Certain combinations of alleles in the Yellowbelly series produce similar phenotypes (e.g.spark vs specter and gravel vs asphalt).Crosses involving these alleles therefore generate progeny whose genotypes are uncertain.Breeders typically label such animals according to their best guess of the animal's genotype.This practice leads to mislabeling by owners and likely explains the apparent mix-ups of animals in our study: one spark/yellowbelly animal was mislabeled by its owner as specter/yellowbelly and vice versa (Fig. 4d); two gravel/yellowbelly animals were mislabeled by their owners as asphalt/yellowbelly (Fig. 5d); and one spark heterozygote was mislabeled by its owner as a yellowbelly heterozygote (Fig. 3c).An additional cause of mislabeling is that breeders are incentivized to label animals as heterozygous for alleles in the Yellowbelly series, even when the phenotype of heterozygotes is difficult to recognize (Fig. 1e  and f).This incentive likely explains why approximately one in four animals described as yellowbelly heterozygotes carried only putatively wildtype copies of EDNRB2 (Fig. 3c).We propose that mislabeling can be reduced in the future by genotyping animals for the EDNRB2 variants reported in our study.

Evolution of endothelin signaling and stripe formation
The phenotypes associated with EDNRB2 variants in ball pythons suggest that endothelin signaling may be required for the development of melanophores, xanthophores, and possibly iridophores in this species.Endothelin signaling in snakes may therefore play a role similar to endothelin signaling in axolotl, where the endothelin ligand EDN3 is thought to direct the development of all three types of chromatophores (Woodcock et al. 2017).Similar regulation has been described for mammals and birds, where endothelin signals regulate the development of melanophores, which are the only type of chromatophore present in these taxa (Saldana-Caboverde and Kos 2010; Bondurand et al. 2018).These species contrast with zebrafish, where endothelin signaling directly promotes the development of iridophores, but is dispensable for xanthophores and plays only an indirect role for melanophores (Frohnhöfer et al. 2013;Patterson and Parichy 2013;Patterson et al. 2014).These observations have led to the proposal that endothelin signaling in vertebrates was ancestrally required for all three types of chromatophores and became specialized for iridophores in teleost fish (Spiewak et al. 2018).Our findings support this hypothesis by adding ball pythons to the list of species in which endothelin signaling likely regulates the development of multiple types of chromatophores.
Dorsal, longitudinal stripes are not unique to morphs in the Yellowbelly series.Other striped morphs in ball pythons include the Genetic Stripe morph, the Champagne morph, and some morphs in the Blue-Eyed Leucistic series (McCurley 2005;Broghammer 2019).Outside ball pythons, stripes occur naturally in many snake species, typically in species that escape predation via rapid flight (Allen et al. 2013).Stripes have arisen multiple times within the snake lineage (Wolf and Werner 1994), and a few species are naturally polymorphic for striped and nonstriped forms (Zweifel 1981;Brodie 1992;King 2003;Westphal and Morgan 2010;Kuriyama et al. 2013;Murakami et al. 2014).Whether stripe formation in these cases involves changes to endothelin signaling or chromatophore number remains to be determined.Our study suggests that dorsal, longitudinal stripes may be easily accessible to the snake body form via genetic changes that alter simple developmental parameters in the neural crest, such as cell number.Candidate genes for controlling these parameters include other genes in the endothelin pathway (Braasch et al. 2009) and genes responsible for specifying cell fates in the neural crest (Mort et al. 2015;Parichy and Spiewak 2015).

Fig. 1 .
Fig. 1.Phenotypes of color morphs in the Yellowbelly series.a) Schematic of homozygote phenotypes in the Yellowbelly series.b) Ball pythons described as homozygous for each allele in the Yellowbelly series.c) Wildtype ball python.Image is used with permission from Brown et al. (2022).d) Ball pythons described as compound heterozygous for alleles in the Yellowbelly series.e) Underbelly of a wildtype ball python.f) Ball python described as heterozygous for the yellowbelly allele.The underbelly is subtly lighter and differently patterned compared to wildtype.b-f) Images are representative examples.Variation exists within each genotype.Photos, Chiron Graves, Gavin Costello, Ball Python Shed, Ball In Hand Pythons, Barbara Penyak, Freedom Breeder, Patrick Faulkner, Off The Wall Balls, Colorado Cold Bloods, Laura Carter, Brent McKelvey, Travis Wyman.

Fig. 2 .
Fig. 2. Phenotypes in the Yellowbelly series are hypothesized to arise from loss of chromatophores.a) Simplified model of chromatophore migration and stripe formation during the larval-to-adult transition in zebrafish (Singh and Nuesslein-Volhard 2015; Patterson and Parichy 2019).Xanthophore precursors persist in the epidermis from larvalhood.Melanophore and iridophore precursors migrate from the neural crest to the epidermis along nerve tracks.The primary route for iridophores is through the horizontal midline.The adult pattern consists of dark stripes (melanophores and sparse iridophores) and light interstripes (xanthophores and dense iridophores).b) Schematic of mutant zebrafish with a primary loss of iridophores and a secondary loss of melanophores(Parichy et al. 2000;Lopes et al. 2008;Frohnhöfer et al. 2013;Patterson and Parichy 2013;Mo et al. 2017;Spiewak et al. 2018).Iridophores fail to populate the skin and remain clustered near the horizontal midline.Melanophores require iridophores for survival and persist only near the iridophores.c) Hypothesized model for chromatophore migration in ball pythons.Chromatophores are presumed to migrate in a dorsal-to-ventral path, as they do in other higher vertebrates(Mort et al. 2015).d-f)Phenotypes in the Yellowbelly series are hypothesized to be caused by a loss of chromatophores.Loss is hypothesized to range from mild (altered patterning) to moderate (striped) to severe (white).This model depicts a loss of all three types of chromatophores, but the loss of iridophores is uncertain.The range of developmental processes plausibly affected includes chromatophore specification, proliferation, migration, or a combination thereof.

Fig. 3 .
Fig. 3.The yellowbelly allele is associated with frameshift variants in EDNRB2.a) Ball python described as homozygous for the yellowbelly allele.Photo, Barbara Penyak.b) Schematic of molecular EDNRB2 alleles associated with the yellowbelly allele.One allele contains a 1-bp deletion.The other allele contains a 17-bp deletion.c) Tableshowingassociation between the yellowbelly allele and the 1-bp and 17-bp deletions.This table includes the original yellowbelly homozygote in which the 1-bp deletion was identified.Numbers of owners and US states indicate the total numbers of owners and US states from which the animals were collected.*, presumed mislabeling by owner, including one animal heterozygous for L152F (see Discussion).d) Schematic of wildtype EDNRB2 protein, created using Protter(Omasits et al. 2014).e, f) Schematics of truncated EDNRB2 proteins resulting from the 1-bp and 17-bp deletions.Schematics are intended to show premature terminations and do not imply that the truncated proteins fold properly and are inserted into the cell membrane.

Fig. 4 .
Fig. 4. The spark and specter alleles are associated with missense variants in EDNRB2.a,b) Ball pythons described as homozygous for the spark allele or the specter allele.Photos, Freedom Breeder, Patrick Faulkner.c) Schematic of molecular EDNRB2 alleles associated with the spark and specter alleles.The spark allele contains missense variant L152F.The specter allele contains missense variant R315P.d) Tableshowingassociations between the spark and specter alleles and the missense variants.This table includes the original spark/yellowbelly compound heterozygote and the original specter homozygote in which the missense variants were identified.Numbers of owners and US states indicate the total numbers of owners and US states from which the animals were collected.*, presumed mislabeling by owner (see Discussion).e, f) Schematics of EDNRB2 proteins containing missense variants, created using Protter(Omasits et al. 2014).g) Alignment of EDNRB2 and EDNRA protein sequences surrounding missense variants.EDNRA is included to demonstrate conservation outside EDNRB2.Full EDNRA sequences were not available for ball python or painted turtle.

Fig. 5 .
Fig. 5.The gravel and asphalt alleles are associated with nucleotide substitutions in EDNRB2 splice donors.a,b) Ball pythons described as homozygous for the gravel allele or the asphalt allele.Photos, Off the Wall Balls, Ball Python Shed.c) Schematic of molecular EDNRB2 alleles associated with the gravel and asphalt alleles.The gravel and asphalt alleles contain G-to-A nucleotide substitution in the −1 positions of the first and sixth splice donors, respectively (SD1 and SD6).d) Tableshowingassociations between the gravel and asphalt alleles and the nucleotide substitutions in splice donors.This table includes the original gravel and asphalt homozygotes in which the nucleotide substitutions were identified.Numbers of owners and US states indicate the total numbers of owners and US states from which the animals were collected.*, presumed mislabeling by owner (see Discussion).e) Alignment of nucleotide sequences for the first and sixth splice donors of EDNRB2.

Fig. 6 .
Fig.6.Summary of genetic variants in EDNRB2 associated with phenotypes in the Yellowbelly series.The most severe phenotype in the Yellowbelly series was all-white skin, which was associated with frameshift variants (1-bp and 17-bp deletions).The next most severe phenotype was dorsal striping, which was associated with missense variants L152F and R315P.The mildest phenotype was altered patterning, which was associated with G-to-A nucleotide substitutions in the −1 positions of splice donors.