how snakes lost their legs?

Report

Loss and Re-emergence o

f Legs in Snakes by Modular Evolution of Sonic hedgehog and HOXD Enhancers

Graphical Abstract

Highlights

d Python legs are truncated due to early arrest of Sonic

hedgehog (SHH) transcription

d The python SHH limb enhancer has weak activity due to

deletion of key binding sites

d HOXD digit enhancers and the footplate expression domain

are conserved in pythons

d Python leg buds form transitory condensations of the tibia,

fibula, and footplate

Leal & Cohn, 2016, Current Biology 26, 2966–2973 November 7, 2016 ª 2016 Elsevier Ltd. http://dx.doi.org/10.1016/j.cub.2016.09.020

Authors

Francisca Leal, Martin J. Cohn

Correspondence [email protected]

In Brief

Leal and Cohn show that hindlimb

development arrests in pythons due to

mutations in an enhancer that controls

Sonic hedgehog transcription in limb

buds. In contrast, HOXD limb enhancers

and distal expression are conserved, and

pythons form a transitory foot skeleton,

providing insights into how legs were lost

and then regained in snake evolution.

Accession Numbers

KX778812

KX778839

KX824111

KX824112

Current Biology

Report

Loss and Re-emergence of Legs in Snakes by Modular Evolution of Sonic hedgehog and HOXD Enhancers Francisca Leal1,2 and Martin J. Cohn1,2,3,4,* 1Howard Hughes Medical Institute 2Department of Biology 3Department of Molecular Genetics and Microbiology UF Genetics Institute, University of Florida, P.O. Box 103610, University of Florida, Gainesville, FL 32610, USA 4Lead Contact

*Correspondence: [email protected]

http://dx.doi.org/10.1016/j.cub.2016.09.020

SUMMARY

Limb reduction and loss are hallmarks of snake evo- lution. Although advanced snakes are completely limbless, basal and intermediate snakes retain pelvic girdles and small rudiments of the femur. Moreover, legs may have re-emerged in extinct snake lineages [1–5], suggesting that the mechanisms of limb devel- opment were not completely lost in snakes. Here we report that hindlimb development arrests in py- thon embryos as a result of mutations that abolish essential transcription factor binding sites in the limb-specific enhancer of Sonic hedgehog (SHH). Consequently, SHH transcription is weak and tran- sient in python hindlimb buds, leading to early termi- nation of a genetic circuit that drives limb outgrowth. Our results suggest that degenerate evolution of the SHH limb enhancer played a role in reduction of hindlimbs during snake evolution. By contrast, HOXD digit enhancers are conserved in pythons, and HOXD gene expression in the hindlimb buds pro- gresses to the distal phase, forming an autopodial (digit) domain. Python hindlimb buds then develop transitory pre-chondrogenic condensations of the tibia, fibula, and footplate, raising the possibility that re-emergence of hindlimbs during snake evolution did not require de novo re-evolution of lost structures but instead could have resulted from persistence of embryonic legs.

RESULTS AND DISCUSSION

Disruption of the SHH/GREM1/AER-FGF Circuit by Early Arrest of SHH Transcription in Python Leg Buds During snake evolution, pythons and boas diverged from the

lineage leading to advanced snakes before hindlimbs were

completely eliminated [6–8]. Python embryos initiate formation

of hindlimb buds, but leg development is not sustained, result-

ing in formation of a rudimentary femur and terminal claw (Fig-

2966 Current Biology 26, 2966–2973, November 7, 2016 ª 2016 Else

ure 1A). We showed previously that python hindlimb buds lack

two critical signaling regions: the zone of polarizing activity

(ZPA) and the apical ectodermal ridge (AER) [9]. The ZPA

controls outgrowth and anteroposterior patterning of vertebrate

limbs by secretion of Sonic hedgehog (SHH) protein [10, 11].

SHH protein was not detected in python hindlimb buds

1–2 days after oviposition [9], and others reported that the

cis-regulatory element that directs limb-specific expression of

SHH was lost in snakes [12, 13]. However, SHH can be acti-

vated in python hindlimb bud cells transplanted under the

AER of chick wing buds [9], suggesting that the mechanism

that controls SHH expression in limbs was not completely

lost in pythons. To resolve this paradox, we first asked whether

SHH transcription occurs in python hindlimb buds at stages

earlier than those examined previously. In Python regius em-

bryos harvested before oviposition and at stage 1 (the day of

oviposition), SHH mRNA was detected in a small domain of

cells at the posterior margin of the hindlimb buds, but this

ZPA-like expression pattern is transient, disappearing within

24 hr of oviposition (Figures 1C and S1A). Analysis of hindlimb

buds in comparably staged anole lizards (Anolis sagrei), which

develop pentadactyl hindlimbs (Figure 1B), revealed a strong

posterior domain of shh that persists after SHH becomes un-

detectable in python hindlimb buds (Figures 1C and 1D). To

test whether weak and transient SHH expression in python hin-

dlimbs results in hedgehog signal transduction, we examined

expression of PTCH1 and GLI1, two transcriptional readouts

of SHH signaling [14]. PTCH1 and GLI1 are expressed in stage

1 python hindlimb buds (Figures 1E and 1G), but the expression

domains are smaller and weaker than those observed in anoles

(Figures 1F and 1H), and expression fades after termination of

SHH expression (Figures 1E and 1G). In limbed tetrapods, GLI3

is expressed anteriorly and distally in limb buds, where it regu-

lates anteroposterior patterning of the digits by repressing SHH

[15]. GLI3 patterns are initially similar in python (stage 1) and

anole (stage 5) hindlimb buds, showing strong anterior-distal

expression that fades near the SHH domain, although at later

stages, GLI3 expression extends further proximally in pythons

(Figures 1I and 1J). Thus, in python hindlimb buds, SHH tran-

scription is initiated in a small group of posterior mesenchymal

cells, and signal transduction occurs; however, SHH expres-

sion is transient, disappearing within 24 hr of oviposition.

vier Ltd.

Figure 1. The SHH/GREM1/AER-FGF

Circuit Is Activated but Not Maintained in

Python Hindlimb Buds

(A and B) Optical projection tomography scans

showing hindlimb skeletal anatomy in a ball python

hatchling (A) and a stage 13 green anole lizard (B).

The axial skeleton is shown in gray, pelvic girdle in

blue, and hindlimb skeleton in red.

(C–L) Gene expression during hindlimb develop-

ment in stage-matched python (C, E, G, I, and K)

and anole (D, F, H, J, and L) embryos at three

stages of development. Broken lines in (I) and (J)

mark proximal limits of GLI3 domains.

(M and N) Apoptosis in python (M) and anole (N)

hindlimb buds stained with LysoTracker Red.

Black arrows in (C) indicate orientation of limb

axes. An, anterior; Po, posterior; Pr, proximal; Di,

distal. See also Figures S1 and S2.

Cessation of SHH transcription is followed by loss of target

gene expression in the posterior region of python hindlimb

buds.

SHH expression in vertebrate limbs is regulated in part by the

AER, which secretes fibroblast growth factors (FGFs), and the

dorsal ectoderm, which produces WNT7a [16–18]. We tested

whether diminished activity of SHH in the python ZPA could

result from deficiencies in either of these two ectodermal

signaling regions. Analysis of WNT7a revealed dorsally compart-

mentalized expression in hindlimb ectoderm (Figure S1G), as

occurs in other limbed tetrapods [16]. Together with our finding

Current Biolog

that LMX1 and EN1 are dorsoventrally

restricted in python hindlimbs [9], this in-

dicates the presence of a dorsal ecto-

dermal signaling region and dorsoventral

polarity in python hindlimb buds. In

limbed tetrapods, AER cells along the

distal edge of the limb buds undergo

pseudostratification and produce FGFs,

which maintain expression of SHH [17,

18]. SHH, in turn, feeds back to maintain

expression of FGFs in the AER, establish-

ing a positive feedback loop that coordi-

nates limb outgrowth and patterning

[19]. A morphological AER is evident in

python hindlimb buds at stage 1 (Figures

S2A and S2B), and FGF8 is expressed in

the AER before oviposition and at stage 1

(Figures 1K and S1B). Shortly after termi-

nation of SHH, the FGF8 domain begins

to degrade, disappearing from the poste-

rior AER between stages 2 and 3, when

weak and patchy expression can be de-

tected anteriorly (Figure 1K). The poste-

rior-to-anterior loss of FGF8 in python

hindlimb buds coincides with the poste-

rior-to-anterior flattening of the AER (Fig-

ures S2C –S2E). The changes in python

AER structure and FGF8 expression

resemble those seen in mouse hindlimb

buds after early deletion of Shh [11], although the severity of

ridge degeneration is greater in pythons, possibly reflecting

disruption of factors in addition to SHH. Anole hindlimb buds,

by contrast, maintain a pronounced AER that expresses Fgf8

throughout (Figure 1L).

In limbed tetrapods, SHH maintains the AER by inducing

Gremlin1 (GREM1), which counteracts the inhibitory activity of

bone morphogenetic proteins (BMPs) on the AER [20]. GREM1

is expressed throughout the anteroposterior axis of the python

hindlimb bud at stage 1 but then weakens following loss of

SHH (Figure S1I). By contrast, expression of BMP4 and its target

y 26, 2966–2973, November 7, 2016 2967

Figure 2. Analysis of trans- and cis-Regula-

tors of SHH Reveals Degeneration and Hy-

pofunctionalization of the Python ZRS

(A–F) Expression of HAND2 (A and B), HOXD13

(C and D), and HOXA13 (E and F) in stage-

matched hindlimb buds of python (A, C, and E) and

anole (B, D, and F) embryos.

(G and H) VISTA conservation plots using human

reference sequence to compare architecture

of the ZRS locus. Peaks indicate conservation >

50%, and colored peaks indicate conservation >

75%. (G) Comparison of ZRS locus conserva-

tion in python, limbed reptiles, and mouse. (H)

Intra-squamate comparison of ZRS locus con-

servation in python, boa, three completely limb-

less advanced snakes, and anole lizard.

(I and J) Functional analysis of preZRS-ZRS en-

hancers from anole (I) and python (J) in transgenic

mice. Columns show, from left to right, whole

embryos at E11.5, forelimb buds, hindlimb buds,

and high-magnification views of the ZPA in fore-

limb (FL) and hindlimb (HL) buds. (I) Anole preZRS-

ZRS drove strong LacZ expression in ZPA of

forelimbs and hindlimbs. (J) Three examples of

transgenic mouse embryos with the python

preZRS-ZRS construct showing weak reporter

activity in forelimbs and hindlimbs. See Experi-

mental Procedures for sample sizes.

(K) Quantification of python and anole PreZRS-

ZRS activity in a luciferase reporter assay dem-

onstrates that the python enhancers have reduced

transcriptional regulatory activity (asterisk in-

dicates significant difference in two-tailed t test,

n = 4, p = 1.6 3 10�4). Error bars indicate standard deviations.

See also Figure S4.

gene, MSX2, is sustained in the distal mesenchyme from stages

1 to 3 (Figures S1J and S1K). In chick limbs, FGF signaling from

the AER drives distal limb development, in part, by regulating cell

survival [21, 22]. Comparison of apoptotic patterns in hindlimb

buds of stage-matched python and anole embryos revealed a

posterior-distal domain of apoptosis in pythons that was not

observed in anoles (Figures 1M and 1N). As the python AER de-

generates during stages 2 and 3, the apoptotic domain expands

distally (Figure 1M), consistent with loss of AER signaling activity.

Together, these results indicate that after the transient pulse of

SHH expression at stage 1, the SHH/GREM1/AER-FGF feed-

back loop breaks down in python hindlimb buds.

Degenerate Evolution of the Python ZRS Underlies Diminished Transcription of SHH Transcription of SHH in tetrapod limbs is activated by binding of

transcription factors to a cis-regulatory element known as the

ZPA Regulatory Sequence (ZRS), a limb-specific enhancer

2968 Current Biology 26, 2966–2973, November 7, 2016

located �1 Mb upstream of SHH, in intron 5 of the LMBR1 gene [23–26]. We first

asked whether the arrest of SHH

transcription in python hindlimbs is asso-

ciated with disruption of trans-acting reg-

ulators of the ZRS. HAND2, HOXD13, and

HOXA13, which encode proteins that

bind the ZRS in mice [23–25, 27], showed strikingly similar

expression patterns in early hindlimb buds of python (pre-ovipo-

sition and stage 1) and anole (stage 5) embryos (Figures 2A–2F,

S1C, and S1D). After the loss of SHH expression in python hin-

dlimb buds, HAND2 expression became progressively weaker

(stages 2 and 3 in Figure 2A), consistent with the finding that

SHH maintains Hand2 posteriorly in mice [28, 29]. HOXA13

and HOXD13 expression persists distally in python (stages 2

and 3) and anole (stages 6 and 7) hindlimbs (Figures 2C–2F).

Although species-specific differences were observed after the

loss of SHH in pythons, HAND2, HOXD13, and HOXA13 expres-

sion were highly conserved in early python and anole hindlimb

buds.

We next asked whether diminished SHH expression in python

hindlimb buds could reflect changes to the integrity of the python

ZRS. Cloning of the �9 kb intron 5 of LMBR1 showed that, in contrast to previous reports that the ZRS has been lost in snakes

[13], pythons have a conserved region that corresponds to the

ZRS of limbed tetrapods (Figure 2G). We also identified a

conserved preZRS, an additional SHH limb enhancer near the

ZRS [30] (Figure 2G). Comparative genomic analysis showed

that despite an overall high degree of similarity in the ZRS and

preZRS sequences in pythons and limbed amniotes, the 50 end of python ZRS contains a region of divergence (Figure 2G),

raising the possibility that mutations in the ZRS could play a

role in diminished activity of SHH in pythons.

We then compared preZRS-ZRS sequences from anole,

python, boa (Boa constrictor, a boid), pit viper (Prothobothrops

mucrosquamatus, a viperid), garter snake (Thamnophis sirtalis,

a colubrid), and king cobra (Ophiophagus hannah, an elapid).

The results show a pattern of degenerative evolution in which

advanced snakes (pit viper, garter snake, and king cobra), which

are completely limbless, show significantly less conservation of

the preZRS and the ZRS than either python or boa, which retain

limb rudiments (Figure 2H). Python and boa have nearly identical

patterns of sequence conservation in the preZRS and the ZRS,

with each showing limited divergence at the 50 end of the ZRS (Figure 2H). Pit viper and garter snake showed less conservation

of the preZRS and the ZRS when compared to python and boa.

The pattern of preZRS conservation was similar in pit viper and

garter snake, but garter snake ZRS showed marked divergence

(Figure 2H). King cobra showed the most extreme degradation of

both enhancers; there was no signature of a preZRS and only

limited conservation at the 30 end of the ZRS (Figure 2H). In order to determine the effects of python ZRS sequence

divergence on its regulatory activity, we generated LacZ reporter

constructs containing the python preZRS-ZRS (preZRS-ZRS-

LacZ) and examined their activity in transgenic mice. A second

reporter construct containing the anole preZRS-ZRS was gener-

ated as a squamate control. The anole preZRS-ZRS drives tran-

scription in a strong ZPA-like domain in mouse forelimbs and

hindlimbs at embryonic day 11.5 (Figure 2I). In contrast, the py-

thon preZRS-ZRS showed minimal activity in mouse forelimbs

and hindlimbs, with patterns ranging from a small number of

LacZ-positive cells to no detectable activity (Figure 2J). To quan-

tify differences in the regulatory activity of the python and anole

preZRS-ZRS, we cloned each enhancer into a luciferase reporter

vector and transfected them into mouse fibroblasts (Figure 2K).

Python preZRS-ZRS showed a 40% reduction in luciferase

transcriptional activation compared to the anole preZRS-ZRS

(Figure 2K). Thus, in vivo and in vitro analyses show that python

preZRS-ZRS is a weak driver of transcription, raising the possi-

bility that divergence of the 50 end of this enhancer underlies diminished SHH expression in python hindlimb buds.

The ZRS is divisible into two domains with different regulato-

ry roles. The 50 end carries spatiotemporal (ZPA-specific) regu- latory information, whereas the 30 end regulates long-range interaction with the SHH promoter [31]. Transcriptional regula-

tors of SHH in the mouse ZPA, such as HAND2 and posterior

HOXD proteins, bind to the 50 end of the ZRS domain [23–25]. We characterized the divergent sequence at the 50

end of the python ZRS and identified three major deletions

(DA, DB, and DC) not found in limbed amniotes (Figure 3A).

The three deletions are shared by three python species

(P. regius, P. molurus, and P. reticulatus) and a boid (Boa

constrictor), indicating their presence in the last common

ancestor of Pythonidae and Boidae.

To determine how each mutation in the python ZRS affects its

activity, we used site-directed mutagenesis to introduce the

mutations individually (DA, DB, or DC) and together (DABC)

into the mouse preZRS-ZRS sequence and then cloned these

into a luciferase reporter construct to assay their activities in

cell culture (Figure S3A). Comparison of the activity of the python

ZRS deletion constructs and a WT mouse ZRS control construct

showed that each of the three python deletions caused a reduc-

tion of luciferase activity, with DC causing the largest decrease

(19.6% of control activity level), followed by DA (53.2%) and

DB (78.8%) (Figure 3B). The ZRS bearing all three python

deletions (DABC) decreased luciferase activity to 12.8% of the

WT control level, a level lower than any of the three individual

deletions (Figure 3B). If the activity that we observed in vitro is

representative of the endogenous effects of these deletions,

then the results suggest that each of the three mutations could

have reduced the efficiency of the ZRS during snake evolution

and that these effects were likely cumulative.

Mutations in Python ZRS Disrupt HOXD Binding Sites Required for Enhancer Activation We next investigated the mechanism by which the DA, DB, and

DC mutations cause reductions in python ZRS activity. In mouse

limbs, the ZRS is transactivated by HAND2, ETS1, and HOXD

proteins [23–25, 32]. HAND2 and ETS1 binding sites have been

characterized at the nucleotide level in mice, and analysis of

these sites in the python ZRS revealed that the HAND2 binding

site is not disrupted but that one of the five ETS1 binding sites

was eliminated (Figure S4). Because the precise positions at

which HOXD proteins bind the 50 end of the mouse ZRS are less well understood, it was unclear how the python ZRS muta-

tions affect its transactivation by HOXD9, HOXD10, or HOXD13.

We addressed this by co-transfecting mouse fibroblasts with a

WT mouse preZRS-ZRS construct or a mouse preZRS-ZRS

construct bearing the python ZRS mutations individually (DA,

DB, or DC) or together (DABC), along with Hoxd9, Hoxd10,

or Hoxd13 expression vectors. The WT mouse preZRS-ZRS

construct showed greatest transactivation by HOXD13 (7.6-

fold), followed by HOXD9 (3.7-fold) and then HOXD10 (1.2-fold)

(Figures 3C and S3B). When python ZRS deletion constructs

with mutations DA, DB, DC, or DABC were co-transfected with

Hoxd9, Hoxd10, or Hoxd13, each of the python deletion con-

structs showed reduced transactivation relative to the mouse

control (Figures 3C, 3D, and S3B). Transactivation was weakest

in the ZRS bearing all three python mutations. For example,

when co-transfected with Hoxd9, Hoxd10, and Hoxd13

expression vectors, the python DABC-ZRS construct showed

only 2% of the activity of the control mouse ZRS that lacked

these mutations (Figures 3C). Thus, each mutation in the python

ZRS weakens its response to HOXD proteins, and all three mu-

tations virtually abolish HOXD transactivation of the python ZRS.

The results described above suggested that HOXD binding

sites were disrupted during evolution of the python ZRS. To

test this hypothesis directly, we asked whether HOXD13,

HOXD10, and HOXD9 can bind to the �400 bp region at the 50 end of the mouse ZRS that corresponds to the domain containing

the three python deletions. Electrophoretic mobility shift assays

(EMSAs) show that all three HOXD proteins bind to this region

of the mouse ZRS (Figure S3C). To determine whether python

Current Biology 26, 2966–2973, November 7, 2016 2969

Figure 3. Molecular Evolution of the Python

ZRS: Loss of Three Binding Sites Abolishes

Transactivation by HOXD Proteins

(A) The 50 end of the python ZRS shows three specific deletions, DA, DB, and DC (DC also has an

8 nt microduplication). These mutations are

conserved in three python and one boid species

(also see Figure S4).

(B–D) Luciferase reporter analyses of mouse

preZRS-ZRS constructs harboring python de-

letions individually (DA, DB, or DC) and together

(DABC). Asterisks indicate significant differences,

two-tailed t test, n = 4. Error bars indicate standard

deviations. See Figure S3B for numerical data and

p values. (B) Effects of python deletions on ZRS

activity. (C) Effects of python mutations on ZRS

transactivation by HOXD13, HOXD10, and

HOXD9. (D) Effects of python mutations on ZRS

transactivation by different combinations of

HOXD13, HOXD10, and HOXD9.

(E) HOXD13-3xFlag binds WT control ZRS regions

that correspond to the deleted regions in the

python ZRS (see Supplemental Experimental

Procedures for oligonucleotide sequences). Each

oligonucleotide corresponds to WT mouse

sequence at positions equivalent to each python

deletion (DA in region A, DB in region B, DC in

region C).

NE, nuclear extracts; NS, non-specific band; black

arrowheads, shifts; red arrowheads, supershifts.

See also Figures S3 and S4.

ZRS mutations DA, DB, or DC reside in a region necessary for

HOXD binding, we designed three short (40 nt) oligonucleotides

from regions of the mouse ZRS that were deleted by mutations

A, B, and C in pythons (Figure S4). Only HOXD13 binds all three

2970 Current Biology 26, 2966–2973, November 7, 2016

oligonucleotides (Figure 3E), suggesting

that python ZRS mutations A, B, and C

each occurred in HOXD13 binding sites

and that its synergy with HOXD9 and

HOXD10 could be mediated by protein-

protein interactions rather than by direct

binding of HOXD10 or HOXD9 to these

regions of the ZRS. These findings,

together with results from our functional

studies of each python ZRS mutation,

show that python mutations DA, DB, and

DC each removed a HOXD13 binding

site that is essential for transactivation of

the ZRS. Thus, by abolishing sites

required for binding of HOXD13, the

mutations at the 50 end of the ZRS can account for the diminished transcription

of SHH in python hindlimb buds.

Conservation of the Distal HOXD Regulatory Landscape and Autopodial Development in Pythons Despite the premature breakdown of the

SHH/GREM1/AER-FGF feedback loop

in python hindlimb buds, our analysis of HOX gene expression

suggested the onset of a late/distal phase of HOXD13 expres-

sion at stage 3 (Figure 2C). In light of the relationship between

distal expression of HOX13 paralogs and digit development in

limbed tetrapods, we monitored the progression of HOXD13 and

HOXA13 in python hindlimb buds at later stages. Surprisingly,

HOXD13 expression spread throughout the distal region of the

python hindlimb bud, overlapping with HOXA13 expression in

a pattern that resembles the autopodial (digit-forming) domains

of HOXD13 and HOXA13 in limbed tetrapods (Figures 4A and

4B) [33, 34]. This late/distal phase of HOXD13 expression was

unexpected given the absence of any remnant of a foot in

pythons.

In mice, distal limb expression of Hoxd13 is controlled by a

series of enhancers in a Topologically Associating Domain

(TAD) centromeric to the HOXD cluster [35]. Analysis of the

syntenic region in python revealed striking conservation of the

distal limb and genital enhancer sequences (Prox, CsA, CsB,

and I-V regulatory islands) that regulate HOXD expression

during tetrapod digit and external genital development (Fig-

ure 4F). The presence of distal HOXD enhancers and the auto-

pod-like patterns of HOXD13 and HOXA13 expression in python

hindlimb buds prompted us to look for evidence of distal hindlimb

skeletogenesis. Analysis of SOX9 expression, which marks

pre-chondrogenic skeletal condensations, revealed a Y-shaped

domain in the proximal and middle region of the hindlimb and a

distal domain that overlaps with the distal region of HOXD13

and HOXA13 expression at stages 4 and 5/6 (Figures 4A–4C).

Over the next two stages, SOX9 expression in python hindlimbs

delineates discrete skeletal condensations that resemble the

three major segments (zeugopod, stylopod, and autopod) of

the tetrapod limb (Figure 4D). Comparison of the SOX9

domains in python in anole hindlimbs suggests that the python

condensations could be anlagen of the femur, tibia and

fibula, and digital plate (Figures 4D and 4E). The fate of

these zeugopodial and autopodial skeletal condensations is

unknown, but they are transitory structures, as distal hindlimb

elements are not found in mature pythons. Taken together, our

results show that pythons have accumulated degenerative muta-

tions in the ZRS that cause precocious arrest of SHH expression

in hindlimb buds, whereas HOXD digit- and genital-specific en-

hancers have been maintained in pythons, and this likely under-

lies the distal activation of HOXD13 and perhaps the specification

of distal skeletal elements in python hindlimb buds.

Conclusions The results presented here suggest that diminished expres-

sion of SHH in the early limb buds of python embryos is a

consequence of three mutations (DA, DB, and DC) in the 50

end of the ZRS that disrupt sites required for ZRS transactiva-

tion by HOXD proteins. We note that our analysis also uncov-

ered a deletion of a single ETS1 binding site in the ZRS, and

although loss of a single ETS1 site is not sufficient to alter

ZRS activity in mice [32], we cannot exclude a role for this

deletion in the diminished activity of the python ZRS. Our

identification of the same three mutations in the ZRS of

pythons and boas indicates that DA, DB, and DC were present

in their last common ancestor (Figure 4G). Based on recent

calibrations of the python and boid clades [6–8], all three

mutations likely arose by the late Upper Cretaceous, when

snakes underwent a major adaptive radiation [1–5, 7, 8]. We

propose that divergence of the ZRS sequence during snake

evolution initially rendered the enhancer hypofunctional,

compromising its ability to drive SHH transcription and, ulti-

mately, resulting in cessation of hindlimb outgrowth and loss

of distal structures (Figure 4G). Furthermore, our observation

that mutations DA, DB, and DC have cumulative effects on

ZRS activity suggests that increased sequence divergence

could have resulted in progressive reduction of the hindlimb

skeleton in basal and intermediate snakes. Hindlimbs disap-

peared completely in advanced snakes (caenophidians), and

comparison of preZRS-ZRS structure across snakes revealed

a pattern of degenerative evolution affecting the preZRS and

the ZRS in vipers, colubrids, and elapids. Furthermore, only

the 30 end of the ZRS is conserved among all snake taxa examined here, suggesting that ZRS degradation began at

the 50 end. Conservation at the 30 end of the ZRS in limbless snakes raises questions about the significance of this region

outside the context of limb development.

Despite the accumulation of mutations in the ZRS and the

reduced nature of python hindlimbs, the genomic and tran-

scriptional machinery necessary to develop limbs has been

largely conserved, from the genetic circuitry active in early

limb buds to the specification of an autopodial (toe-forming)

domain. Conservation of HOXD enhancers centromeric to the

HOXD cluster in snakes suggests that this regulatory domain

remained under selection after its role in digit development

became obsolete. Our recent finding that HOXD genes, but

not SHH, are transcribed in the developing hemipenes of py-

thons [36], together with evidence that HOXD expression in

digits and external genitalia are under shared regulatory control

[37], suggests that the HOXD distal enhancers were retained in

snakes due to their essential role in external genital develop-

ment. Finally, our discovery that python embryos develop

transitory cartilage condensations of the lower leg and foot

suggests that re-acquisition of fully developed hindlimbs in

extinct snakes [1–5] may not have required de novo re-evolu-

tion of lost structures, but could have resulted from persistence

of the embryonic legs.

EXPERIMENTAL PROCEDURES

Animal protocols were reviewed and approved by the University of Florida

Institutional Animal Care and Use Committee. Full experimental procedures

can be found in the Supplemental Experimental Procedures.

ACCESSION NUMBERS

The accession numbers for the sequences reported in this paper are GenBank:

KX778812-KX778839 and KX824111-KX824112.

SUPPLEMENTAL INFORMATION

Supplemental Information includes Supplemental Experimental Procedures

and four figures and can be found with this article online at http://dx.doi.org/

10.1016/j.cub.2016.09.020.

A video abstract is available at http://dx.doi.org/10.1016/j.cub.2016.09.

020#mmc3.

AUTHOR CONTRIBUTIONS

F.L. and M.J.C. designed the experiments. F.L. performed the experiments.

F.L. and M.J.C. analyzed the data, interpreted the results, and wrote the

manuscript.

Current Biology 26, 2966–2973, November 7, 2016 2971

Figure 4. ZRS Degeneration and Cryptic

Hindlimb Development in Python Embryos

Suggests a Model for Loss and Re-acquisi-

tion of Legs during Snake Evolution

(A and B) Python hindlimb buds exhibit a late distal

autopodial phase of HOXD13 and HOXA13

expression.

(C–E) SOX9 expression delineates prechondro-

genic condensations in three discrete regions

(stylopod, zeugopod, and autopod) of stage-

matched python (C and D) and anole (E) hindlimbs.

Note that SOX9 distal domain develops within the

autopodial domain of HOXD13/HOXA13 expres-

sion in python (compare A and B with C and D).

fm, femur; fb fibula; ad, autopodial domain. Black

arrows show limb axes: An, anterior; Pr, proximal;

Di, distal; Po, posterior.

(F) VISTA conservation plots (using human

sequence as reference) show that python and boa

retain the conserved regulatory archipelago of

HOXD digit and genital enhancers in the gene

desert centromeric to the HOXD cluster. Peaks

indicate the presence of conserved sequences A

(CsA) and CsB of the Global Control Region, Prox,

and the five regulatory islands.

(G) Model for the role of degenerate evolution of

the ZRS in snake hindlimb reduction. Phylogeny

adapted from [3, 5, 6, 8]. Hindlimb structures (red)

shown at right were generated by optical projec-

tion tomography (except for fossil species). ZRS is

indicated by green (functional), yellow (hypofunc-

tional), or red (degenerated); solid outlines indicate

known ZRS status; broken outlines indicate pre-

dicted status. Tree branches are black for lineages

with known ZRS structure, gray for predicted

structure, and broken for extinct lineages. (1)

Ancestral ZRS (green) drove SHH expression and

digit development in limbed squamates. (2) We

propose that after the divergence of Serpentes

from limbed squamates, mutations DA, DB, and

DC weakened ZRS function (yellow), reducing

SHH activity and causing loss of distal hindlimb

elements. Yellow/green ZRS reflects uncertainty

about whether digit reduction occurred in the

fossil Najash or whether absence of digits reflects

taphonomic conditions. (3) Absence of a tibia,

fibula, foot, and toes in scolecophidians (e.g., blind

snakes) and alethinophidians suggests ZRS

hypofunctionalization in the common ancestor.

Transitory prechondrogenic condensations of

zeugopodial and autopodial skeleton form in em-

bryos but then degenerate, resulting in absence of

these structures in adults. (4) Re-emergence of the

tibia and fibula (zeugopod) and digits (autopod) in

Tethyan snakes may have resulted from per-

sistence of embryonic skeletal condensations

into adulthood. The requirement of the ZRS for

digit development suggests that ZRS activity may

have been amplified, despite the presence of

deletions. (5) Further degeneration of the ZRS (red) in advanced snakes (caenophidians) rendered it non-functional (see Figure 2H), resulting in total limblessness.

Skeletal morphologies pictured next to the tree terminals were derived from OPT scans, except for fossil species.

ACKNOWLEDGMENTS

We are grateful to K. Kelley and K. Backer-Kelley at the UF ICBR electron

microscopy core for assistance with SEM; J. Wisby, B. Cole, O. Tarazona,

B. Caruso, T. Sanger, and B. Kircher for assistance with animal husbandry

and/or eggs; B. Cole and Colección Herpetológica Universidad

2972 Current Biology 26, 2966–2973, November 7, 2016

Industrial de Santander (UISR-2549) for adult snake material; R. Renne

and H.-S. Choi for help with luciferase assays; and C. Larkins and

D. Menke for sharing materials. F.L. thanks the Gans Collections and

Charitable Fund. This project was supported the Howard Hughes

Medical Institute (M.J.C). F.L. is an HHMI international student research

fellow.

Received: July 14, 2016

Revised: August 29, 2016

Accepted: September 12, 2016

Published: October 20, 2016

REFERENCES

1. Caldwell, M.W., and Lee, M.S.Y. (1997). A snake with legs from the marine

Cretaceous of the Middle East. Nature 386, 705–709.

2. Rage, J.C., and Escuillié, F. (2000). Un nouveau serpent bipède du

Cénomanien (Crétacé). Implications phylétiques. C. R. Acad. Sci. Paris

Earth Sci. 330, 513–520.

3. Houssaye, A., Xu, F., Helfen, L., Buffrénil, V.D., Baumbach, T., and

Tafforeau, P. (2011). Three-Dimensional Pelvis and Limb Anatomy of the

Cenomanian Hind-Limbed Snake Eupodophis descouensi (Squamata,

Ophidia) Revealed by Synchrotron-Radiation Computed Laminography.

J. Vertebr. Paleontol. 31, 2–7.

4. Tchernov, E., Rieppel, O., Zaher, H., Polcyn, M.J., and Jacobs, L.L. (2000).

A fossil snake with limbs. Science 287, 2010–2012.

5. Apesteguı́a, S., and Zaher, H. (2006). A Cretaceous terrestrial snake with

robust hindlimbs and a sacrum. Nature 440, 1037–1040.

6. Pyron, R.A., Burbrink, F.T., and Wiens, J.J. (2013). A phylogeny and

revised classification of Squamata, including 4161 species of lizards and

snakes. BMC Evol. Biol. 13, 93.

7. Head, J.J. (2015). Fossil calibration dates for molecular phylogenetic

analysis of snakes 1: Serpentes, Alethinophidia, Boidae, Pythonidae.

Palaeontol. Electronica 18.1.6FC, 1–17.

8. Zheng, Y., and Wiens, J.J. (2016). Combining phylogenomic and superma-

trix approaches, and a time-calibrated phylogeny for squamate reptiles

(lizards and snakes) based on 52 genes and 4162 species. Mol.

Phylogenet. Evol. 94 (Pt B), 537–547.

9. Cohn, M.J., and Tickle, C. (1999). Developmental basis of limblessness

and axial patterning in snakes. Nature 399, 474–479.

10. Riddle, R.D., Johnson, R.L., Laufer, E., and Tabin, C. (1993). Sonic hedge-

hog mediates the polarizing activity of the ZPA. Cell 75, 1401–1416.

11. Zhu, J., Nakamura, E., Nguyen, M.-T., Bao, X., Akiyama, H., and Mackem,

S. (2008). Uncoupling Sonic hedgehog control of pattern and expansion of

the developing limb bud. Dev. Cell 14, 624–632.

12. Sagai, T., Hosoya, M., Mizushina, Y., Tamura, M., and Shiroishi, T. (2005).

Elimination of a long-range cis-regulatory module causes complete loss of

limb-specific Shh expression and truncation of the mouse limb.

Development 132, 797–803.

13. Sagai, T., Masuya, H., Tamura, M., Shimizu, K., Yada, Y., Wakana, S.,

Gondo, Y., Noda, T., and Shiroishi, T. (2004). Phylogenetic conservation

of a limb-specific, cis-acting regulator of Sonic hedgehog ( Shh).

Mamm. Genome 15, 23–34.

14. Pearse, R.V., 2nd, Vogan, K.J., and Tabin, C.J. (2001). Ptc1 and Ptc2 tran-

scripts provide distinct readouts of Hedgehog signaling activity during

chick embryogenesis. Dev. Biol. 239, 15–29.

15. Büscher, D., Bosse, B., Heymer, J., and Rüther, U. (1997). Evidence for

genetic control of Sonic hedgehog by Gli3 in mouse limb development.

Mech. Dev. 62, 175–182.

16. Yang, Y., and Niswander, L. (1995). Interaction between the signaling

molecules WNT7a and SHH during vertebrate limb development: dorsal

signals regulate anteroposterior patterning. Cell 80, 939–947.

17. Laufer, E., Nelson, C.E., Johnson, R.L., Morgan, B.A., and Tabin, C. (1994).

Sonic hedgehog and Fgf-4 act through a signaling cascade and feedback

loop to integrate growth and patterning of the developing limb bud. Cell 79,

993–1003.

18. Niswander, L., Tickle, C., Vogel, A., Booth, I., and Martin, G.R. (1993).

FGF-4 replaces the apical ectodermal ridge and directs outgrowth and

patterning of the limb. Cell 75, 579–587.

19. Niswander, L., Jeffrey, S., Martin, G.R., and Tickle, C. (1994). A positive

feedback loop coordinates growth and patterning in the vertebrate limb.

Nature 371, 609–612.

20. Scherz, P.J., Harfe, B.D., McMahon, A.P., and Tabin, C.J. (2004). The limb

bud Shh-Fgf feedback loop is terminated by expansion of former ZPA

cells. Science 305, 396–399.

21. Rowe, D.A., Cairns, J.M., and Fallon, J.F. (1982). Spatial and temporal

patterns of cell death in limb bud mesoderm after apical ectodermal ridge

removal. Dev. Biol. 93, 83–91.

22. Fallon, J.F., López, A., Ros, M.A., Savage, M.P., Olwin, B.B., and Simandl,

B.K. (1994). FGF-2: apical ectodermal ridge growth signal for chick limb

development. Science 264, 104–107.

23. Galli, A., Robay, D., Osterwalder, M., Bao, X., Bénazet, J.D., Tariq, M.,

Paro, R., Mackem, S., and Zeller, R. (2010). Distinct roles of Hand2 in initi-

ating polarity and posterior Shh expression during the onset of mouse limb

bud development. PLoS Genet. 6, e1000901.

24. Osterwalder, M., Speziale, D., Shoukry, M., Mohan, R., Ivanek, R., Kohler,

M., Beisel, C., Wen, X., Scales, S.J., Christoffels, V.M., et al. (2014).

HAND2 targets define a network of transcriptional regulators that

compartmentalize the early limb bud mesenchyme. Dev. Cell 31, 345–357.

25. Capellini, T.D., Di Giacomo, G., Salsi, V., Brendolan, A., Ferretti, E.,

Srivastava,D., Zappavigna,V., and Selleri,L.(2006). Pbx1/Pbx2 requirement

for distal limb patterning is mediated by the hierarchical control of Hox gene

spatial distribution and Shh expression. Development 133, 2263–2273.

26. Lettice, L.A., Heaney, S.J., Purdie, L.A., Li, L., de Beer, P., Oostra,

B.A., Goode, D., Elgar, G., Hill, R.E., and de Graaff, E. (2003).

A long-range Shh enhancer regulates expression in the developing

limb and fin and is associated with preaxial polydactyly. Hum. Mol.

Genet. 12, 1725–1735.

27. Johnson, E.J., Neely, D.M., Dunn, I.C., and Davey, M.G. (2014). Direct

functional consequences of ZRS enhancer mutation combine with

secondary long range SHH signalling effects to cause preaxial poly-

dactyly. Dev. Biol. 392, 209–220.

28. Charité, J., McFadden, D.G., and Olson, E.N. (2000). The bHLH tran-

scription factor dHAND controls Sonic hedgehog expression and estab-

lishment of the zone of polarizing activity during limb development.

Development 127, 2461–2470.

29. Litingtung, Y., Dahn, R.D., Li, Y., Fallon, J.F., and Chiang, C. (2002). Shh

and Gli3 are dispensable for limb skeleton formation but regulate digit

number and identity. Nature 418, 979–983.

30. Park, K., Kang, J., Subedi, K.P., Ha, J.H., and Park, C. (2008). Canine poly-

dactyl mutations with heterogeneous origin in the conserved intronic

sequence of LMBR1. Genetics 179, 2163–2172.

31. Lettice, L.A., Williamson, I., Devenney, P.S., Kilanowski, F., Dorin, J., and

Hill, R.E. (2014). Development of five digits is controlled by a bipartite long-

range cis-regulator. Development 141, 1715–1725.

32. Lettice, L.A., Williamson, I., Wiltshire, J.H., Peluso, S., Devenney, P.S., Hill,

A.E., Essafi, A., Hagman, J., Mort, R., Grimes, G., et al. (2012). Opposing

functions of the ETS factor family define Shh spatial expression in limb

buds and underlie polydactyly. Dev. Cell 22, 459–467.

33. Nelson, C.E., Morgan, B.A., Burke, A.C., Laufer, E., DiMambro, E.,

Murtaugh, L.C., Gonzales, E., Tessarollo, L., Parada, L.F., and Tabin, C.

(1996). Analysis of Hox gene expression in the chick limb bud.

Development 122, 1449–1466.

34. Zakany, J., and Duboule, D. (2007). The role of Hox genes during verte-

brate limb development. Curr. Opin. Genet. Dev. 17, 359–366.

35. Montavon, T., Soshnikova, N., Mascrez, B., Joye, E., Thevenet, L.,

Splinter, E., de Laat, W., Spitz, F., and Duboule, D. (2011). A regulatory ar-

chipelago controls Hox genes transcription in digits. Cell 147, 1132–1145.

36. Leal, F., and Cohn, M.J. (2015). Development of hemipenes in the ball

python snake Python regius. Sex Dev. 9, 6–20.

37. Lonfat, N., Montavon, T., Darbellay, F., Gitto, S., and Duboule, D. (2014).

Convergent evolution of complex regulatory landscapes and pleiotropy

at Hox loci. Science 346, 1004–1006.

Current Biology 26, 2966–2973, November 7, 2016 2973