Introduction

Evolutionary transitions have always fascinated researchers as they pose important questions about the developmental and genetic origins of evolutionarily novel traits (Shubin et al. 2009). Due to (relatively) easy access to relevant fossils and the availability of tractable experimental models, one of the best studied evolutionary transitions is the origin and evolution of the tetrapod limb (Shubin et al. 1997; Shubin 2002; Schneider and Shubin 2013).

A century of descriptive, experimental and theoretical research has produced unprecedented insights into the evolution of the tetrapod limbs. The emergence and adoption of molecular techniques have made it possible to start to reveal the identity of key genes in the process and to discern the dynamics of the gene regulatory networks that underpin morphological changes in the vertebrate appendage. Building on these emerging data, multiple (broadly similar) hypotheses have been put forward to explain how the fins of ancestral tetrapods have gradually changed and adapted to walk on land. Decades of collaborative research have resulted in an evolutionary framework that explains, in broad brushstrokes, how molecular changes can lead to the appearance of long bones and the evolution of the autopod (Capdevila and Belmonte 2000; Wagner and Chiu 2001; Freitas et al. 2014; Tanaka 2016; Amaral and Schneider 2018; Lalonde and Akimenko 2018), and we can even try to answer questions such as why do we have five digits (Tabin 1992; Wagner and Chiu 2001; Saxena et al. 2017).

Furthermore, innovative and imaginative experiments have also provided important insights into the evolution of phenotypic variability in extant tetrapod groups, from limblessness in reptiles to the evolution of powered flight in bats (Cooper et al. 2012; Leal and Cohn 2018; Royle et al. 2021).

In this brief review, we summarize our knowledge about the development of the stereotypical tetrapod limb and present the current experimental evo-devo framework used to analyze its evolution. Finally, we also argue that the easy access to exponentially growing genomic resources and widespread use of recently developed molecular techniques that can be applied to virtually any species of interest herald a new era of evo-devo research.

Overview of tetrapod limb development

From the early pioneering works of Ross G. Harrison and Viktor Hamburger through those of John Saunders and Lewis Wolpert, the past century has revolutionized our understanding of the tetrapod limb development (Oppenheimer 1966; Kirk and Allen 2001; Tickle 2017; Smith 2021). As several recent high-profile reviews have summarized the essentials about the emerging “textbook view” of this process (Tickle 2015; Petit et al. 2017; McQueen and Towers 2020; Royle et al. 2021), we will only briefly introduce the basic anatomical structure of the tetrapod limbs and, subsequently, the main molecular players involved (see Table 1).

Table 1 A summary of the most important genes involved in vertebrate limb development
Full size table

The limbs of tetrapods are partitioned into three major endochondral segments (Fig. 1A). The stylopod (humerus in the forelimb and femur in the hindlimb), found on the most proximal side, is an elongated tubular bone, articulating with the pelvic girdle. A more distal segment, the zeugopod, consists of two longer tubular bones (radius and ulna in the forelimbs, tibia and fibula in the hindlimbs). The most distal and exclusive portion of tetrapod appendages is the autopod, a short tubular bony structure that involves the wrist or ankle, respectively, and the fingers or toes. This basic common design of “one bone [that] articulates with two bones, which attach to a series of small blobs, which connect with the fingers or toes” underlies the architecture of all tetrapod limbs (Shubin 2008).

Fig. 1
figure 1

Major phases in the development of the tetrapod appendage. A The schematic structure of a stereotypical tetrapod limb. Limbs consist of three major segments: the proximal stylopod (purple), the zeugopod (orange and green), and the distal-most autopod (red, cyan and yellow). B The specification of limb identity during embryogenesis. The location of both the fore-, and hindlimb buds is determined by a rostrocaudal “Hox code”. Rostral Hox proteins directly activate the expression of Tbx5, thus promote forelimb identity. Caudal Hox factors support the formation of hindlimbs through Pitx1 and Tbx4 expression. C Phases of embryonic limb development with respect to the interactions between the main regulators. In the initiation stage, the Fgf-secreting AER cells promote the emergence of the posteriorly located ZPA, which expresses Shh. The morphogen Shh activates the expression of the BMP-antagonist Grem1 in the mesenchymal tissue of the developing limb bud. Grem1 will ultimately maintain Fgf expression through BMP suppression. During propagation, the AER-derived Fgf-level increases, which induces the suppression of Grem1. As the result of this negative feedback, BMP levels will be stabilized in the forming limb, thus the Fgf expression is eventually silenced, which leads to the termination of embryonic limb development. D The role of 5’ HoxD gene paralogs during limb development. The early phase of HoxD expression (labeled with green) regulated by 3’ located enhancers will organize the development of proximal segments of the limb (i.e., stylopod and zeugopod). In this stage, the Hoxd8-d11 paralogs are the most active. The second or late phase of HoxD expression (marked with yellow) is essential for autopod formation. The expression of Hoxd10-d13 paralogs is controlled by 5’ situated regulator elements centered on the posterior side of the limb bud. E Digit formation is orchestrated by a Turing-like mechanism. The striping pattern of chondrogenic digit precursors is created by interactions between Sox9-, Wnt-, and BMP-signaling. The expression pattern of Sox9, the master regulator of cartilage development overlaps with the location of future digits and is regulated by the aforementioned signaling pathways, AER-derived Fgf signals and distally expressed Hoxd13. The closer the Sox9 expression is located to the AER, the higher the wavelength of the expression pattern becomes. (Based on Raspopovic et al., 2014). F Formation of the DV axis in the developing limbs. The mesenchymal Lmx1b expression is regulated by ectodermal Wnt7a, and it is essential for the formation of the dorsal side. On the ventral side of the developing limb bud, BMP signals induce the expression of ectodermal En1 and to the procuration of ventral identity

Full size image

Limb development is initiated by the localized activation of fibroblast growth factor (Fgf) signaling pathway, more precisely by expression of the Fgf10 gene encoding one of the pathway’s ligands in the lateral plate mesoderm (LPM) of the developing body wall (Jin et al. 2019; Royle et al. 2021). Through the induction of canonical Wnt-signaling (in particular, the expression of Wnt3a) in the surface ectoderm, Fgf-signaling establishes the apical ectodermal ridge (AER) (Kengaku et al. 1998), a signaling center that will be essential for the development of the limb bud (Fig. 1B, C). The AER itself is a source of Fgf signals (its cells secrete Fgf8 (Mahmood et al. 1995)), and its activity will contribute to the proliferation of underlying cells, promoting the outgrowth of the limb bud (McQueen and Towers 2020).

Current models suggest that the position of limb initiation is determined by a “Hox code” along the anteroposterior (AP) axis. Homeobox (Hox) genes can indirectly regulate the expression of T-box transcription factor 4 (Tbx4) and T-box transcription factor 5 (Tbx5), master regulators of limb development and identity, with Tbx5 being responsible for forelimb development and the concerted action of Tbx4 and Paired-like homeodomain 1 (Pitx1) defining hindlimb identities (Fig. 1B) (Logan and Tabin 1999; Takeuchi et al. 1999; Royle et al. 2021).

As the limb grows, the AER will promote the emergence of a second signaling center, the zone of polarizing activity (ZPA), in the posterior part of the bud. The ZPA secretes the morphogen Sonic hedgehog (Shh), which not only regulates the growth of the bud, but also provides the positional information essential for establishing the AP axis of the bud (Tickle et al. 1975; Riddle et al. 1993).

ZPA-derived Shh will indirectly antagonize Bone morphogenetic protein (BMP) signaling in this region. This is achieved by the induction the expression of the BMP-antagonist Gremlin 1 (Grem1) in the limb mesenchyme (Khokha et al. 2003). Grem1 will contribute to sustained Fgf-signaling, as BMPs would suppress Fgf-signaling. During the growth of the limb, however, a negative feedback between high Fgf levels and Grem1 expression will ultimately induce and stabilize the expression of BMPs in the limb (Pignatti et al. 2014). Consequently, Fgf signals from the AER will be silenced, and in their absence the ZPA also disappears, terminating the embryonic phase of limb development (Fig. 1C) (Tickle and Towers 2017).

The role of Hox genes in limb development is not limited to limb initiation, though. For example, genes of the HoxD cluster are expressed in multiple waves in the limb bud, regulated by distinct regulatory elements, and ultimately these genes will play an essential role first by establishing stylopod, zeugopod and autopod identities, and later by contributing to digit development (Zákány et al. 2004).

Concomitant mutation of paralogous 5’ genes in the HoxA and HoxD clusters have revealed the roles of these genes in establishing proximo-distal (PD) identities in the developing limb. Double mutants for Hox10, Hox11 and Hox13 paralogs will display abnormalities in the size and shape of their stylopods, zeugopods and autopods, respectively (Davis et al. 1995; Fromental-Ramain et al. 1996; Wellik and Capecchi 2003).

Pioneering work from Denis Duboule’s group has demonstrated that this early expression of the HoxD cluster is regulated by 3’ situated early regulatory elements, whereas the later, autopod-specific expression of the Hoxd10-13 genes will depend on the activity of 5’situated global control region (GCR) (Zákány et al. 2004; Andrey et al. 2013). The GCRs of the autopod progenitors are regulated by the Shh morphogen gradient emanating from the ZPA, with high Shh concentrations inducing the expression of all the aforementioned Hoxd genes, and progressively lower morphogen concentrations turning off Hoxd10, -11 and -12 expression, respectively, so that the future thumb expresses only Hoxd13 (Fig. 1D) (Zákány et al. 2004).

For a long time, the development of the fingers has been one of the more mysterious aspects of tetrapod limb development. Lately, however, through the elegant combination of experimental and modeling work the Sharpe group have clearly demonstrated that a Turing-like mechanism is responsible for creating the chondrogenic pattern of finger development (Sheth et al. 2012; Raspopovic et al. 2014).

Alan Turing’s original “reaction–diffusion” model posited that the expression of a particular diffusible substance that can promote its own expression as well as the expression of its more diffusible repressor will create characteristic patterns from a homogeneous state (Turing 1952).

Earlier experimental work has already established that in the chondrogenesis of the finger bones multiple signaling pathways play a role, but were not able to reveal neither how their interplay will lead to the stereotypical “five-digit” pattern of the tetrapod autopod, nor how the mutations of individual components will result in the observed malformations. The model from the Sharpe group suggests that the main molecular players of this developmental step create an intricate, three-component “reaction–diffusion” system, which is also differentially modulated along the PD axis. According to the model, a combination of Hoxd13, expressed in the limb mesenchyme, and Fgf signals, emanating from the AER will modulate the interplay between Wnt- and BMP-signaling and the expression of Sox9, the master regulator of cartilage development (Fig. 1E). The combination of inductive events with negative feedback loops and auto-inhibition will lead to the pattern of mesenchymal condensations that is observed during finger development (Raspopovic et al. 2014).

Besides the formation of the PD and AP axis, development of the limb also coincides with the establishment of the dorso-ventral (DV) identities (Fig. 1F). In the developing bud, the dorsal mesenchymal expression of LIM homeobox transcription factor 1 beta (Lmx1b) is essential to establish dorsal identities, as Lmx1b mutant mice will grow footpads (of ventral identity) on the dorsal part of their paws (Riddle et al. 1995; Vogel et al. 1995). The expression of Lmx1b is regulated by the expression of Wnt7a in the dorsal ectoderm (Parr and McMahon 1995). In contrast, in the ventral side BMP signals will lead to the induction of Engrailed 1 (En1) in the ectoderm, which in turn will suppress Wnt7a expression. Accordingly, the impairment of BMP-signaling in the limb bud will lead to the loss of En1 expression and the appearance of dorsal identities on both sides of the limb (Pizette et al. 2001; McQueen and Towers 2020).

Gene expression changes during the fin-to-limb transition

While both the development of fins and those of tetrapod limbs is orchestrated by a very conserved set of genes, the later stages of development are strikingly distinct in these two appendage types. Digits, the most distal elements of the tetrapod autopod, were long considered as an evolutionary novelty, not homologous to the dermal fin rays of fish. More recent studies, however, suggest that the development of both dermal fin rays and digits are governed by the same mechanisms regardless of their clear histological differences (Nakamura et al. 2016; Onimaru et al. 2016).

When studying the radiation of sarcopterygians, a conspicuous tendency for the acquisition of novel distal endochondral skeletal elements and the reduction of dermal fin rays can be observed (Fig. 2A). The most important fossils that help to reconstruct the emergence of tetrapod limbs are dated from the second half of the Devonian period. Of the currently known tetrapodomorphs (tetrapod-like sarcopterygians) from this period, the Tiktaalik roseae already had several tetrapod-like characteristics, but instead of fingers it still had long fibrous fin rays (i.e., lepidotrichia) on the most distal region of their appendages (Shubin et al. 2006). A detailed analysis of the recently discovered complete Elpistostege watsoni fossil also revealed that this species belonged to one of the most advanced stem tetrapod groups known. At least two of the multiple digits identified in the forelimb of these animals are composed of a number of aligned, subequal unbranched radials that articulate in a one-to-one relationship. This peculiar structure further blurs the boundary between sarcopterygians and tetrapods (Cloutier et al. 2020). Extant non-tetrapod sarcopterygian species, such as lungfish also possess more distal endoskeletal elements and smaller apical fin folds (AF)—a homolog of the tetrapod AER—compared with fishes that diverged earlier during the evolutionary history (Hodgkinson et al. 2009).

Fig. 2
figure 2

The evolutionary origin of the tetrapod limb based on fossils and molecular data. A Phylogenetic tree of gnathostomes showing the skeletal homologies between different taxa. Well established homologies are color-coded (see Fig. 1 for color legends). Catshark (Scyliorhinus canicula) is a member of the basal gnathostome group, the cartilaginous fishes. Bichir (Polyodon spathula) shows basal ray-finned fish characteristics (i.e., all three segments can be found in the pectoral fins), whereas the more derived teleosts, such as zebrafish (Danio rerio), possess solely anterior and distal radials. Sarcopterygians are represented by eight different taxa, including two extant fish species (Latimeria, Neocaratodus forsteri), one extinct fish (Eusthenopteron), two tetrapodomorphs (Tiktaalik and Elpistostege), and one basal and two modern tetrapods (Acanthostega, Tulepreton and Mus musculus, respectively). Elpistostege had at least two digits, which are homologous to tetrapod ones. The basal Acanthostega showed archaic polydactylous limbs, whereas more derived tetrapods possesses maximum five fingers. Note that the lepidotrichia of Elpistostege is not included in this figure. B Comparison between the early stages of appendage development of teleost fish and tetrapods. During fin bud development, the apical ectodermal ridge undergoes folding, creating the AF. The cells of the AF will give rise to connective tissues of the fin, including the long dermal fin rays. Due to the folding of the apical ectodermal ridge mesenchymal cells proximal to the AF receive less Fgf signals, they proliferate less and the endochondral elements of the fin become reduced. The late phase of hoxd13 expression is also restricted to more distal regions of the bud. In contrast, the AER of the tetrapod limb buds will provide ample proliferative signals for the nearby mesenchymal cells and will induce the robust expression of Hoxd13 and other 5’ Hox genes. The proliferating mesenchymal cells will ultimately give rise to developing long and structurally complex endochondral bones

Full size image

An important hallmark of the transition has been the remodeling of the fin rays into distal endoskeletal elements that most likely required the stepwise acquisition of a number of exaptations, co-options and loss-of-function mutations in genes related to fish fin development. Indeed, some genes [e.g., actinodin 1 (and1) and actinodin 2 (and2)] that encode structural proteins of fin rays are completely absent from tetrapod genomes (Zhang et al. 2010). During fin development, the apical ectodermal cells fold, creating the AF, thus able to give rise to connective tissues of the fin blade. However, in tetrapod appendages, the AER cells do not undergo folding and therefore mesenchymal cells are capable of proliferating as a response to AER signals, creating endochondral skeletal elements of the autopod (Dudley et al. 2002).

Loss of and1 and and2 function in zebrafish larvae results in reduced fins and less developed (or absent) AF in the case of pectoral fins (Zhang et al. 2010). If we accept that the reduction of AF size had significant consequences for vertebrate appendage morphology and might have been a prerequisite for later limb evolution, it is tempting to consider that the loss of these two genes was very important during the fin-to-limb transition.

Posterior HoxA and HoxD genes are essential for normal development of appendages in both fish and tetrapods. In the developing fins of the zebrafish (Danio rerio), the expression of posterior hox genes shows the same temporal expression pattern as observed in early tetrapod limbs (Ahn and Ho 2008). Late-phase Hox expression governs the growth of digits in tetrapod and dermal fin rays as well as in the distal radials in fish (Nakamura et al. 2016; Onimaru et al. 2016). Furthermore, triple mutant zebrafish for hoxa13a, hoxa13b and hoxd13a have shortened dermal fin rays and possess supernumerary distal radials (Nakamura et al. 2016).

On the other hand, induced overexpression of hoxd13a in 30–32 h post-fertilization (hpf) zebrafish larvae (a stage just prior to the anterior expansion of late-phase 5’hox gene expression in fins) results in the reduction of the AF (Freitas et al. 2012). In the absence of AF, the cells of the distal chondrogenic tissue are able to expand more distally, causing increased cell proliferation, which leads to the formation of an autopod-like outgrowth. In pectoral fins, hoxd13a overexpression also affects the expression patterns of other genes involved in appendage development. For instance, fgf8 and and1 show reduced expression domains, while hoxa13, ETS variant transcription factor 4 (etv4) and cytochrome P450, family 26 (cyp26), all involved in the proximo-distal subdivision of the appendages, are upregulated. Overall, these expression patterns are highly similar to that observed in tetrapod limbs (Freitas et al. 2012).

The cis-regulatory logic of Hox genes appears to be almost identical in both groups, as many of the regulatory sequences are conserved between fish and tetrapods (Gehrke et al. 2015). However, comparison between mice, chicken, zebrafish, and gar, a basal ray-finned fish that lacks the teleost-specific whole-genome duplication (TGD) revealed that some enhancers are peculiar to tetrapods, implying that new sequences must have been acquired for autopod evolution. Interestingly, many cis-elements of HoxA/D genes are conserved from basal fishes to tetrapods, but not in teleosts, which suggests that many enhancers were lost or changed beyond recognition due to the teleost-specific TGD. It is also noteworthy that gar (but not zebrafish) enhancers are active in transgenic mice during development, indicating that a subset of trans-regulatory elements were already present in the last common ancestor of actinopterygians and sarcopterygians (Gehrke et al. 2015).

Formation of limbs with multiple long bones, separated by joints, another hallmark feature of the tetrapod autopod might have required the upregulation of specific signaling pathways. Recently identified gain-of-function mutations in the zebrafish genes wiskott-aldrich syndrome-like b (waslb) and vav guanine nucleotide exchange factor 2 (vav2) both result in the development of supernumerary bones (termed “intermediate” radials). While it is not clear if these supernumerary bones are indeed orthologous to any of the long bones of the tetrapod limb, further observations have demonstrated that Wasl is also required for limb patterning in mice and it can modulate the expression of hoxa11b in the zebrafish fin (Hawkins et al. 2021).

Waslb encodes a regulator of the actin cytoskeleton signaling, while Vav2 encodes a guanine nucleotide exchange factor. The discovery of these genes highlights the fact that our understanding of limb pattern formation is still incomplete and there are molecular pathways yet to be explored that could influence limb development (Hawkins et al. 2021).

Evolution of the stereotypical tetrapod limb

While the emergence of the stereotypical tetrapod appendage with one bone in the stylopod, two bones in the zeugopod and five fingers in the autopod was a real evolutionary milestone, the evolution of limbs did not end there. The radiation of tetrapod groups produced plenty of instances when this aforementioned basic pattern was modified, sometimes even resulting in the loss of the limbs.

It is worth to emphasize that the relative recency of the novel mutations driving the emergence of group-specific limb modifications also means that some of these mutations can be “reverted” as demonstrated by the abundance of limb-related atavistic characters in these groups (Hall 1984).

Next, we will summarize our current knowledge about the gene expression changes related to the evolution of the tetrapod limb.

Limb loss in snakes (and other reptiles)

One of the more extreme changes to the tetrapod body plan can be identified as the loss of limbs observed in snakes. Earlier discoveries have identified the fossil records of transitional forms belonging to this group with four or two appendages (Tchernov et al. 2000; Apesteguía and Zaher 2006; Martill et al. 2015); therefore, we can reconstruct the evolution of this group with reasonable confidence.

Interestingly, the serpentiform body plan has evolved multiple times independently in squamates, suggesting that this group might have several pre-adaptations that facilitate the evolution of limblessness (Leal and Cohn 2018). Evolutionarily the loss of the pectoral girdle and forelimbs preceded the reduction of the hindlimbs as demonstrated not only by the fossil record, but also by the presence of rudimentary pelvic girdle and hindlimb buds in some basal snake lineages (Boughner et al. 2007).

The alterations in gene expression and the resulting developmental changes that have led to the loss of the forelimb are still under investigation. While initial observations suggested that the positional information necessary for the development of pectoral girdle was missing from snakes due changes in the mesenchymal “Hox code” (Cohn and Tickle 1999), later studies revealed that at least in some species, the primaxial expression domains of Hox genes are normal, but Tbx5 shows a deregionalized expression, which might lead to the impaired development of the forelimb (Woltering et al. 2009; Head and Polly 2015).

In stark contrast to the ambiguity surrounding the molecular causes behind missing forelimbs, important advances in recent years have revolutionized our understanding of the reduction of the snake hindlimbs. Almost as soon as the regulatory sequence responsible for the expression of Shh in the developing ZPA has been identified, it was also showed that this highly conserved sequence, present in many vertebrate groups from teleosts to mammals was impaired in snakes (Sagai et al. 2004, 2005). Later research has showed that this enhancer region, the ZPA Regulatory Sequence (ZRS), has accumulated multiple mutations during the evolution of extant snake species and these mutations impair the binding of transactivating Hox and ETS proto-oncogene 1 (Ets1) transcription factors, which leads to hypofunctionality and only weak induction of Shh expression (Kvon et al. 2016; Leal and Cohn 2016). Using an innovative knock-in experimental approach, Kvon and coworkers have also demonstrated that the restoration of a single 17 bp deletion in the hypofunctional python ZRS was sufficient to reinstate the full functionality of this enhancer (Kvon et al. 2016). This important insight also points to a relatively straightforward explanation to the appearance of atavistic hindlimb “revertants” in this group, as it might be relatively easy to restore ZRS functionality (Fig. 3A).

Fig. 3
figure 3

Major changes during the evolution of extant tetrapod limbs. A The role of a limb-specific Shh enhancer, the ZRS (blue), in tetrapods. Transgenic mice carrying snake-derived ZRS showed impaired limb-specific Shh expression, thus reduced limbs compared to wildtype mice. However, restoration of a 17-bp-long region of ZRS in mutant mice successfully rescued the wildtype phenotype, implying a relatively easy explanation of the appearance of revertant snake mutants in nature (based on Kvon et al. 2016). B Two examples for the origin of “webbing” in the interdigital regions in autopods. The developing wings of bats exhibit elevated expression of Grem, which suppresses the Bmp-signaling essential for the apoptosis of cells in the interdigital tissues. The elevated levels of Fgf-signaling in the interdigital regions also impair apoptosis. Somewhat similarly in water fowl, such as ducks, elevated expression of Grem also suppresses Bmp-induced apoptosis. (Based on Weatherbee et al. 2006). C The developmental origins of the reduction of fingers in cows, camels and horses. In cows, a relative decrease in the area of active Shh-signaling results in a reduction in the number of fingers. In contrast, in three-toed jerboas, camels and horses ectopic foci of apoptosis during relatively late stages of limb development result in the elimination of lateral finger primordia. (For autopodal color coding see Fig. 1A) (Based on Cooper et al. 2014 and Lopez-Rios et al. 2014). D Homology between the digits of the stereotypical, five-fingered vertebrate limb (on the left) and the three digits of bird wings (on the right) based on gene expression analysis. Different colors encode different digit identities

Full size image

It is also worth mentioning that analogous changes in the ZRS might also explain the evolution of other limbless reptiles as well (Sagai et al. 2004).

Limb loss and hyperphalangy in cetaceans

Besides reptiles, another prominent group characterized by (hind-)limb loss is the Cetaceans, a group of mammals that forfeited their terrestrial lifestyle and returned to water. This transition is documented in exquisite detail in the fossil record as well (Pyenson 2017); therefore, we have a very good idea about the chronology of the anatomical changes that contributed to the evolution of Cetacean characteristics.

Due to sampling difficulties, it is harder to reconstruct the gene expression changes that contributed to the reduction of the pelvic girdle and the hindlimb. Some pioneering work has identified characteristic changes in the expression of Fgf8 and Shh during the development of the spotted dolphin (Stenella attenuata) and suggested that the hypofunction of these genes is probably behind the transient appearance and later degeneration of the hindlimb bud (Thewissen et al. 2006). The premature downregulation of Fgf8 apparently results in the early loss of the AER, which might also explain why Shh expression cannot be maintained in the ZPA (Thewissen et al. 2006). In the absence of sustained Shh expression, similarly to serpentiform lizards, only a reduced pelvic girdle and hindlimb can develop. This effect might also be helped by the accelerated evolution of Hoxd11 which has acquired some characteristic amino acid substitutions in this lineage (Li et al. 2019).

The evolution of 5’ Hox genes might be also associated with a characteristic feature of the cetacean forelimbs (flippers), namely, hyperphalangy (Cooper et al. 2018). While all other mammalian digits are composed of three phalanges, Cetaceans have evolved supernumerary phalanges in their forelimb to support their aerodynamic flippers. At this point, our understanding about the molecular changes related to this exquisite feature of the Cetacean limb is extremely limited, but it has been noted that both Hoxd12 and Hoxd13 show signs of selection and/or acquired some characteristic mutations in this group (Wang et al. 2009).

In the developing flipper, Fgf8 also appears to be expressed in the interdigital tissues of late fetal stage animals and might contribute for a sustained expression of Gremlin in these cells as well. Based on what we know about the anti-BMP properties of Gremlin in the developing bat wing (see below), this observation might also provide a mechanism for the maintenance of soft interdigital tissues characteristics for cetacean flippers (Cooper et al. 2018).

Gene expression changes during the evolution of the bat wing

Bats (Chiroptera), the only mammalian group capable of powered flight, underwent a series of physiological and morphological adaptations that suit their aerial lifestyle (Cooper et al. 2012). For the purpose of our review, we will focus on the development of wings, which consist of membranous skin stretched between the elongated digits of the forelimb. Interestingly, both aforementioned aspects of the peculiar bat wing morphology have been linked to the modulation of BMP-signaling.

The elongation of digits results from high levels of chondrocyte proliferation during early stages, which can be linked to elevated levels of Bmp2-driven signaling (Sears et al. 2006). The skin membrane stretching between the fingers arises due to defects in the apoptotic process that eliminates the interdigital cells during the development of the limbs in other tetrapods. The inhibition of apoptosis is partly due to persistent elevated expression of Gremlin, a BMP inhibitor that is downregulated during later stages of mouse forelimb development, and as a consequence high BMP levels will trigger cell death (Fig. 3B). Apoptosis is also inhibited by the expression of Fgf8 in the limb mesenchyme, which is an evolutionarily new expression domain as this gene is expressed in the AER under normal circumstances (Weatherbee et al. 2006).

Later research also uncovered some important regulatory changes in the expression of limb-related genes that most likely also contributed to the evolution of the bat wing. For example, analysis of a bat-specific enhancer of the Paired related homeobox 1 (Prx1) homeobox gene, which is important for limb skeletal elongation, also suggests that evolutionarily novel traits can arise due to small changes in multiple regulatory elements (Cretekos et al. 2008).

While this hypothesis seems trivial, in the absence of the adequate high-throughput technologies, earlier studies have been reduced to “candidate-gene” approaches, where only the expression of genes previously implicated in limb development has been analyzed in detail. In the latter years, this has markedly changed and newly developed techniques have been applied to the biological questions related to bat wing evolution as well. A combination of genomic approaches allowed the identification of so-called bat accelerated regions (BARs), which cover rapidly evolving genomic sequences that overlap with predicted limb enhancers of mice (Booker et al. 2016). Of the over 160 BARs identified this way, an enhancer of the HoxD cluster has been also found that might be responsible for the enhanced expression of 5’ Hoxd genes in the developing wing (Booker et al. 2016).

A systematic recent transcriptomic approach has also uncovered a large number of genes that are differentially expressed in the bat forelimb and hindlimb. More importantly, subsets of long non-coding RNAs (lncRNAs) have also been identified that are transcribed as anti-sense transcripts of developmentally important genes (e.g., Tbx5-as1) and show peculiar differences in their expression between the developing bat forelimb and hindlimb (Eckalbar et al. 2016). This unbiased approach also confirmed the upregulation of BMP-antagonists (i.e., Gremlin) and Fgf-signaling components during later stages of forelimb development, in line with the reduction of apoptosis in the interdigital tissues. Notably, it has also been suggested that the downregulation of canonical Wnt-signaling and the upregulation of the Wnt/planar cell polarity (PCP) pathway in the wings might also contribute to the morphological differences between them and the hindlimbs (Eckalbar et al. 2016).

Limb evolution in ungulates and the jerboa

Digit loss has evolved independently in multiple mammalian groups. We can observe it in jerboas, horses and multiple artiodactyl species. Recent elegant work from the Tabin and Zeller labs has uncovered multiple developmental pathways in these species that drive this convergent phenomenon (Fig. 3C) (Cooper et al. 2014; Lopez-Rios et al. 2014).

In some species, such as the camel, the horse and the three-toed jerboa (Dipus sagitta) a reduction of the Fgf8 expression domain in the AER will result in a post-patterning excess of apoptosis in the more proximal-lateral cells of the autopod, eliminating the digits at the edges of the developing limb (Cooper et al. 2014). In contrast, the developing pig and cow embryos show a reduction in the levels of autopodial Shh-signaling, resulting in patterning defects that can explain the reduction in the number of digits (Lopez-Rios et al. 2014).

Jerboas are also noteworthy in the context of limb evolution as they have evolved elongated long bones (tibia and metatarsus) in their hindlimbs, partly explained by enlarged chondrocytes in the growth plates of these bones (Cooper et al. 2013; Moore et al. 2015). Since the size of chondrocytes is regulated by Insulin like growth factor 1 (Igf1) in mice, a modulation in the expression of this gene might also explain the evolution of jerboa hindlimbs; however, this hypothesis has not been tested yet (Cooper et al. 2013).

Evolution of limbs in birds

Despite being one of the central models of experimental developmental biology for a century, chick limb development still has some mysterious aspects. The developing wing bud was essential to elucidate the function of several signaling centers essential for appendage development, yet the seemingly simple question of digit identity in the chicken forelimb has been the subject of intense debate for decades (Young et al. 2011).

Three major hypotheses have been put forward to reconcile anatomical, paleontological and embryological evidences. The Pyramid Reduction Hypothesis (PRH) posits that the identity of digits in the avian wing is equivalent to those of digits II, III and IV in the stereotypical tetrapod limb. The Axis Shift Hypothesis (ASH) argues for a homology between avian digits and tetrapod digits I, II and III. Finally, the Frame Shift Hypothesis (FSH) suggests no straightforward equivalences between the digits of extant birds and their archosaurial ancestors and the identities of digits I, II, II have shifted to a different position in the avian wing (Young et al. 2011).

Recent detailed transcriptomic analyses of bird and reptile developing forelimbs have provided evidence that a combination of PRH and FSH seems to have occurred and identities of the wing digits correspond to that of digits I, III and IV in the plesiomorphic archosaur ancestor (Wang et al. 2011; Stewart et al. 2019).

Digit identity is not the only intriguing evolutionary change that has occurred in birds, though. The loss of flight in ratites, for example, also involved changes in overall wing morphology and development. A recent study related to the wing development in the emu (Dromaius noveahollandie) has found that heterochronic shifts observed in the developing emu wing, compared with those of chicken, can be attributed to a reduction in Fgf10 expression levels and consequently impaired Fgf-signaling, which coincide with a reduction in proliferation (Young et al. 2019). As ectopic expression of Fgf10 in emu wing buds was able to induce precocious limb development, the authors argued that the modulation of Fgf10 expression due to the loss of specific enhancers could be behind the reduction in wing size. However, while open chromatin sequences, specific to the chicken wing bud and absent from the emu could be identified in the vicinity of Fgf10, in an in vivo test this region failed to show enhancer identity. On the other hand, using the same experimental setup a putative enhancer near the Spalt-like transcription factor 1 (Sall1) gene was also identified. Sall-1 is also involved in limb development downstream of Fgf- and Wnt-signaling, and showing high differential expression between emu forelimb and hindlimb could be confirmed as limb-specific enhancer (Young et al. 2019).

The webbing on the feet of waterfowl has also been a relatively well-studied example of evo-devo adaptation. Similarly to later findings in bats, the webbing is related to a lack of apoptosis in the interdigital tissues of the developing hindlimb. In this case, an upregulation in the activity of Gremlin has been linked to the BMP-antagonism necessary to the suppression of cell death (Fig. 3B) (Laufer et al. 1997; Merino et al. 1999; Pizette et al. 2001).

Finally, an analysis of domesticated pigeon breeds (Columba livia) that have feathered legs has showed that the new phenotype can be attributed to a transformation in limb bud identity, as the early hindlimb buds show elevated expression of Tbx5 and reduced expression of Pitx1, thus an intermediate gene expression profile between classical tetrapod fore- and hindlimbs (Domyan et al. 2016). Similar changes in Tbx5 expression could be also detected in the legs of feathered chicken breeds as well.

Methodological advances in the study of limb development and evolution

In the past few decades, our knowledge about the development and evolution of vertebrate appendages has been progressively expanded. A combination of forward and reverse genetic experiments in multiple model organisms have revealed the most important genes involved in these processes.

Genomic approaches have also revealed that the conserved function of multiple vertebrate genes involved in appendage development can be attributed to a large set of conserved enhancers regulating the early stages of embryonic appendage development (Gehrke et al. 2015; Adachi et al. 2016; Lettice et al. 2017).

Our expanding knowledge base has also revealed that while most genes controlling appendage development are highly conserved among all examined gnathostome species, their regulation shows marked differences in some taxa, which leads to a remarkable phenotypic diversity. While traditionally the study of non-coding sequences has been more difficult, in the past two decades the mapping of regulatory sequences has become more feasible due to an explosion in the number of available genome sequences and the emergence of improved in silico and laboratory techniques.

For example, the Vertebrate Genomes Project (VGP) aims to create reference-level sequences for all vertebrate species (Rhie et al. 2021). While this huge undertaking will probably take several years, side projects that aim to sequence the genomes of a huge number of species from taxa relevant for limb evolution as well, such as Bat1K (Teeling et al. 2016), Fish10K (Fan et al. 2020) or Bird10K (Feng et al. 2020) will surely contribute to further important insights into the sequence-level evolution of coding and regulatory sequences.

Modern in silico methods are not only able to locate conserved non-coding regions in genomes, but also capable of predicting transcription factor binding sites within them with high confidential rates (Oshchepkov and Levitsky 2011; Fang et al. 2016; Vandel et al. 2019). Thus, less conserved, taxon- or species-specific enhancers can be mapped, helping to understand the evolutionary background of different appendage morphologies.

For the experimental verification of putative regulatory activities of given non-coding genomic regions several recently developed or improved techniques can be applied. One of the most commonly utilized methods is the Chromosome Conformation Capture (3C) and its derivatives (i.e., 4C, 5C, Hi-C). These methods take advantage of the fact that during the initiation of transcription many cis-regulatory elements physically interact with the promoter of a certain gene, resulting in chromosomal looping. 3C and its derivatives are capable of capturing the frequency of chromosomal interactions within a specific cell type or tissue type in a chosen developmental stage, showing the sequences that of might have regulatory potential (Kempfer and Pombo 2020; McCord et al. 2020).

The mapping of developmentally relevant enhancers could also take advantage of some physical and epigenetic features of these regulatory sequences. During gene activation, the nucleosomal structure of the chromosome becomes relaxed, thus allowing molecules to bind to the promoter and other regulatory sequences. The Assay for Transposase-Accessible Chromatin using sequencing (ATAC-seq) technique is able to map these open regions genome-wide, using a hyperactive bacteria-derived transposase that inserts adapters into accessible loci for subsequent analysis (Buenrostro et al. 2013; Corces et al. 2017).

Chromatin immunoprecipitation with sequencing (ChIP-seq) utilizes antibodies, coupled with a high molecular weight bead for precipitation, specific for certain histone marks, transcription factors, or the RNA Polymerase II (Park 2009). This method is also able to create genome-wide maps of DNA regions marked by specific epigenetic tags, such as acetylation on the lysine 4 residue of histone 3 (H3K4ac), associated with enhancer activity (Calo and Wysocka 2013). One downside of ChIP-seq is that it suffers from low signal-to-noise ratio due to the preceding cross-linking step (Nakato and Sakata 2020). More optimized variants of the technique, such as CUT&RUN and CUT&Tag, bypass the cross-linking, thus are capable of mapping with high signal-to-noise ratios (Skene and Henikoff 2017; Kaya-Okur et al. 2019).

With the combination of the aforementioned techniques, one is able to obtain detailed maps of putative regulatory sequences. However, in order to examine the activity pattern and redundancy of the sequences of interest, more laborious experiments are also required. The spatial and temporal activity of putative enhancers has been traditionally characterized using transposon-mediated (Kawakami 2007; Korzh 2007; Ivics et al. 2009) or site-directed transgenesis techniques (Mosimann et al. 2013; Roberts et al. 2014). The advent of various CRISPR/Cas9 methodologies (Gilbert et al. 2013; Mali et al. 2013; Qi et al. 2013; Lopes et al. 2016) facilitated the emergence of more complex analyses, such as in situ enhancer dissections (Letelier et al. 2018; Hörnblad et al. 2021). Conventional mutagenesis methods usually knock out the entire regulatory region; however, during fine-scale dissection well-defined subregions can be characterized, which is important to understand molecular nature of cis-acting regulation. As elegant studies related to snake ZRS evolution have already demonstrated, the comparative analysis of orthologous enhancers from different taxa using a combination of comparative genomics and genome editing techniques can also help us to discern the evolution of regulatory regions (Kvon et al. 2016).

Conclusions and outlook

Over the past century, the formation of the tetrapod limb became one of the most studied and most important research topics of developmental biology. Key observations made in the context of developing limbs had a profound influence on the main concepts of the field, from the role of morphogen gradients in tissue patterning, to the effect of long-range enhancers on gene function, up to the chromosomal dynamics of Hox clusters (Petit et al. 2017). As non-classical model systems started to gain more popularity in the past few decades and our understanding about the species-specific differences in appendage development exploded, the molecular networks underlying the development of the vertebrate limb also became a central evo-devo paradigm (Gehrke and Shubin 2016; Howenstine et al. 2021).

The recent expansion of the molecular toolbox by the aforementioned high-throughput techniques (i.e., 3C, ATAC-seq, ChIP-seq, etc.), combined with cutting-edge developments, such as single-cell spatial transcriptomics (Srivatsan et al. 2021), will help us add much needed experimental details to current models of limb development. The exponential increase in the available data will also help us refine these theoretical models. Finally, perhaps for the first time in the history of experimental embryology it is feasible to extend these studies to multiple representative species of multiple taxa, thus creating a comprehensive description of past and ongoing events in the evolution of the tetrapod limbs.