- Research article
- Open Access
Transcriptomic insights into the genetic basis of mammalian limb diversity
BMC Evolutionary Biology volume 17, Article number: 86 (2017)
From bat wings to whale flippers, limb diversification has been crucial to the evolutionary success of mammals. We performed the first transcriptome-wide study of limb development in multiple species to explore the hypothesis that mammalian limb diversification has proceeded through the differential expression of conserved shared genes, rather than by major changes to limb patterning. Specifically, we investigated the manner in which the expression of shared genes has evolved within and among mammalian species.
We assembled and compared transcriptomes of bat, mouse, opossum, and pig fore- and hind limbs at the ridge, bud, and paddle stages of development. Results suggest that gene expression patterns exhibit larger variation among species during later than earlier stages of limb development, while within species results are more mixed. Consistent with the former, results also suggest that genes expressed at later developmental stages tend to have a younger evolutionary age than genes expressed at earlier stages. A suite of key limb-patterning genes was identified as being differentially expressed among the homologous limbs of all species. However, only a small subset of shared genes is differentially expressed in the fore- and hind limbs of all examined species. Similarly, a small subset of shared genes is differentially expressed within the fore- and hind limb of a single species and among the forelimbs of different species.
Taken together, results of this study do not support the existence of a phylotypic period of limb development ending at chondrogenesis, but do support the hypothesis that the hierarchical nature of development translates into increasing variation among species as development progresses.
Drivers of morphological diversification include the evolutionary birth of new genes and changes in the regulation of orthologous genes. For some time, many biologists have argued that modifications in gene regulation have triggered many if not most of the evolutionary changes in morphology that characterize the history of life [1–12]. However, the specific changes in gene regulation and expression that have driven morphological diversity remain unknown for most systems. This study employs a comprehensive, multi-species comparison of the gene expression in developing mammalian limbs to investigate how the expression of the highly conserved genes governing organ morphogenesis have been modified during the diversification of limb morphology within and among species.
Mammalian limbs represent an exceptional case study for investigating the evolution of gene expression, as they have undergone extensive evolutionary diversification  and represent a well-characterized model system for development [14, 15]. Furthermore, from the wings of bats to the flippers of whales to the one-toed hooves of horses, mammalian limb diversification has been crucial to the success of the group. For example, through the morphological diversification of their limbs, mammals were able to infiltrate almost every habitat in the world, and exhibit a wide-range of feeding and social behaviors . To date, studies of mammalian limb evolution and development has been limited mostly to investigations of the role of candidate genes (e.g., bats [3, 11, 12, 16–19], whales , opossums [20–23], non-cetacean artiodactyls , jerboas ), or, in a couple of cases, transcriptomes [24–27], in individual species. From these studies we have learned that mammalian limb diversification has proceeded not by major changes to limb structure (e.g., complete loss or gain of entire segments), but by the modification of segments inherited from their generalized, pentadactyl ancestor [28, 29]. Accordingly, these studies of the candidate genes and transcriptomes of isolated species studies suggest that mammalian limb diversification has likely been driven largely by the differential expression of conserved genes shared by all mammals [2–4, 8, 11, 30–32].
While studies of the candidate genes and transcriptomes of isolated species have significantly advanced our understanding of mammalian limb evolution, and morphological evolution in general, there are several outstanding questions concerning limb developmental evolution that have proven difficult to answer with these approaches. For example, we do not yet have a comprehensive view of when, developmentally, the expression pattern of shared genes diverges among the fore- and hind limbs of a single species, or among the limbs of different species. Development is a hierarchical process that builds on itself as time progresses. Because of this, some biologists have proposed that the developmental processes mediating earlier developmental events (e.g., initial specification of organ fields) might be less variable than those governing later events (e.g., organ specialization) [33–40]. The earlier stages of limb development also coincide with the proposed phylotypic period for vertebrates, which is the period in which vertebrate species most closely resemble each other [41, 42]. The phylotypic period has been proposed to encompass the initiation and early development of the limb bud prior to the onset of chondrogenesis [37, 43–45]. Strong inductive signaling between different parts of the embryo has been proposed to characterize the phylotypic period, and result in the conservation of this stage across vertebrates [37, 41, 46]. Evidence for a phylotypic period has recently been observed in other systems, including plants , flies , zebrafish , nematodes , and across several vertebrates . Our first working hypothesis is therefore that the expression pattern of shared genes tends to diverge later (i.e., at or after the onset of chondrogenesis) rather than earlier (i.e., before the onset of chondrogenesis) in limb development.
We also do not have a complete picture of which shared genes and gene pathways differ the most in their expression among the fore- and hind limbs of a single species and among the limbs of different species. Candidate gene approaches have identified differences in the expression of major limb patterning genes (e.g., Shh, Fgf’s, Hox genes, Bmp’s, etc.) across mammalian species [3, 4, 8, 12, 16, 17, 19]. However, these studies were not comprehensive in their gene coverage. Our second working hypothesis is that the expression patterns of several, major limb patterning genes significantly differ between the fore- and hind limbs of a single mammalian species, and among the limbs of different mammalian species.
Finally, we do not know the degree to which different species share a common pattern of gene expression divergence (e.g., timing, genes involved, etc.) between their fore- and hind limbs. Similarly, we do not know the degree of similarity between the patterns of gene expression divergence within (e.g., between the fore- and hind limb of a single species) and among (e.g., between the forelimbs of two different species) species. If certain aspects of the genetic basis of limb development are more evolvable than others, then we might expect the pattern of gene expression divergence to be conserved within species, and within and among species. Our third working hypothesis is therefore that the same genes are differentially expressed in the fore- and hind limbs of several species, and among the limbs of species.
To test these hypotheses we performed the first comprehensive, multi-species comparison of the transcriptomes of developing mammalian limbs. We compared the overall patterns of gene expression in the limbs of four mammals that are taxonomically diverse and represent extremes of mammalian limb development and adult limb structure, namely bats (Carollia perspicillata), pigs (Sus scrofa), mice (Mus musculus), and opossums (Monodelphis domestica). These species include a generalized, pentadactyl mammal (mouse), a representative of the only mammalian group to evolve a wing capable of powered flight (bat), a marsupial whose hind limb development lags significantly behind forelimb development relative to eutherian mammals (opossum) , and a mammal that displays digit reduction (pig). We performed RNA-Seq on the developing fore- and hindlimbs of these four mammalian species . We then quantified the transcriptomes of the fore- and hind limbs of these species at the ridge, bud, and paddle stages of development (Fig. 1) [54–56]. The limb first begins to grow out from the flank during the ridge stage, becomes a semicircle that is as long as it is wide during the bud stage, and forms a hand-plate during the paddle stage. Condensation of the limb cartilaginous elements of the limb begins between the bud and paddle stages [23, 56–58]. We then used the results to determine the pattern of transcriptomic divergence among the forelimbs of different species, the hind limbs of different species, and the fore- and hind limbs of a single species. For several genes, we also assay their spatial expression domains using whole mount in situ hybridization.
Results of the among-species analyses of this study are consistent with the pattern of gene expression during initial limb outgrowth being conserved among mammals relative to subsequent stages of limb development, while within-species results are more mixed. In regard to the former, results further suggest that the development of the limbs of different species likely begins to diverge before the onset of cartilage condensation. Results also suggest that the expression patterns of the limb development genes Hand1, Isl1, Myog, Pax1, Tbx4, Tbx5, and Tnnt2 significantly differ between the fore- and hind limbs of all mammalian species. This study also identified several genes with known roles in limb development that display greatly divergent expression in the fore- and hind limbs of most species and among the limbs of all species (e.g., Col2a1, Hoxa13, Mecom, Pitx1, Rarb, Tbx4, Wnt5a, Zbtb16). However, results suggest that, on the whole, the genes that are differentially expressed in fore- and hind limbs differ from species to species, and that the genes that are that are differentially expressed within (e.g., between the fore- and hind limb of a single species) and among (e.g., between the forelimbs of two different species) species are generally distinct. This study therefore identifies a small subset of genes with possible roles in the generation of the distinct morphologies of the fore- and hind limbs, and a small but different subset with possible roles in the evolutionary divergence of limb morphology across species. This result, combined with the observed differences in the overall patterns of among- and within- species variation, suggests that the processes controlling within- and among- species variation may differ. This study therefore provides tangible insights into the pattern of gene expression divergence during the specialization of mammalian limbs, and the coupled achievement of mammalian diversity.
Timing of the divergence of gene expression of shared genes within species (fore- vs. hind limb) – Limbs at the ridge, bud, and paddle stages (Fig. 1) were removed from each species and total RNA purified (see Materials and Methods and  for more details). We used the Illumina TruSeq RNA Sample Preparation Kit to generate libraries. Once libraries were generated and sequenced, data were quality controlled before analysis (See Additional file 1 and Additional file 2: Figures S1 & Additional file 3: Figure S2 for full details). We generated heat-maps to visualize the similarity of gene expression profiles for the fore- and hind limbs of all developmental stages for a given species (Fig. 2), using the pvclust R package to determine the significance of the clusters  and (Additional file 4: Figure S3). In bats (Fig. 2a), the fore- and hind limbs cluster together at the paddle stage of development, with the forelimb bud falling at the next step out. The forelimb ridge and hind limb bud also cluster together, with the hind limb ridge as the outlier for all limbs and stages. Of these, the clusters of the paddle stage fore- and hind limbs and of the forelimb ridge and hind limb bud are significant (P-value < 0.05) (Additional file 4: Figure S3A). In mouse, the hind limb bud and paddle cluster together with the forelimb paddle as the next step out, and the fore- and hind limb ridges cluster together (Fig. 2b). The mouse forelimb bud is the outlier. The cluster of the fore- and hind limb ridge, fore- and hind limb paddle, and hind limb bud stages is significant (Additional file 3: Figure S2B). In opossum, the fore- and hind limb buds cluster together and the fore- and hind limb ridges cluster together, with the forelimb paddle clustering with the bud samples and the hind limb paddle falling as an outgroup to all other samples (Fig. 2c). None of these clusters are significant (Additional file 4: Figure S3C). Finally, in pig the fore- and hind limbs cluster together at every developmental stage, with the ridge and bud stages clustering together to the exclusion of the paddle samples (Fig. 2d). Of these, the bud and paddle clusters of the fore- and hind limbs are significant (Additional file 4: Figure S3D). Thus the overall results of the clustering are mixed. The ridge staged fore- and hind limbs cluster together in three of the four species, and the bud and paddle staged limbs cluster in two. However, only four of these clusters are statistically supported (one ridge, one bud, and two paddle).
Timing of the divergence of shared gene expression among species (fore- vs. forelimb, hind vs. hind limb) - Pairwise species coefficients were generated for 6,583 orthologs in the 4 species to measure between-species gene expression conservation at each developmental stage (See Additional file 1 for more details). The results of these correlations were used to generate heat-maps to visualize the similarity of gene expression profiles for the forelimbs of all developmental stages for all species, and for the hind limbs (Fig. 3). All species were positively correlated (Spearman coefficient >0.5), however, none of these correlations were significant at the P-value < 0.05 level (Additional file 5: Figure S4). For the ridge stage of both the fore- and hind limb, and the bud stage of the hind limb, bat and pig cluster together, and mouse and opossum cluster together (Fig. 3a, d, e). For pig and bat, this clustering pattern matches the phylogenetic relationships between these species . At the bud stage of forelimb development (Fig. 3b), bat and pig clusters remain clustered together, followed by opossum, and then mouse. Finally, for the paddle stage of development, pig and mouse cluster together, followed by opossum and then bat for the forelimb (Fig. 3c). In the hind limb, mouse and opossum cluster together, followed by pig and then bat (Fig. 3f). Therefore bat, with its highly divergent limb morphology, is an outlier for the paddle stage of development for both the fore- and hind limb.
We also calculated the conservation of the gene expression profiles of bat, mouse, opossum, and pig across embryonic limb development, using the mean of all species pairwise Spearman coefficients. All resulting Spearman coefficients are positive and > 0.50, suggesting that the orthologous genes might perform similar functions between species (Figs. 4a, 4b). In both the fore- and hind limb, the degree of gene expression conservation decreases from the ridge (forelimb = 0.6011, hind limb = 0.6170) to bud (forelimb = 0.5697, hind limb = 0.5882) to paddle (forelimb = 0.5613, hind limb = 0.5799) stage. When all the genes are sampled, the distributions of gene expression conservation levels between the ridge, bud, and paddle stages are significantly different (T-test, P-value < 0.05*). To test the robustness of this difference with respect to the selection of orthologous genes, we randomly sub-sampled 500 sets of orthologous genes at all developmental stages at intensities ranging from 50 to 100% of all orthologous genes. According to the resulting distributions (Fig. 4), only 70% of the genes need to be sampled to find a statistically significant difference (non-overlapping 95% confidence intervals) between the ridge and bud stages of forelimb and of hind limb development. However, the significances of the differences between the bud and paddle stages of the forelimb and of the hind limb are more dependent on the genes that are sampled. The distributions indicate that over 90% of orthologous genes must be sampled to find significant differences (non-overlapping 95% confidence intervals) between these stages.
To further investigate the pattern of gene expression divergence among species, we calculated the evolutionary age of the genes of each species in our study. For this analysis, evolutionary age is defined by the phylogenetic origin of a set of genes, or its phylostratum. Gene sets with lower phylostrata have a younger evolutionary age, and vice versa. We found that bat, pig, mouse, and opossum have 3,828, 10,841, 21,693, and 14,789 genes, respectively, that are specific to a single species (i.e., they do not have a homologue in any other species) (Fig. 5a). For this, all genes (protein-coding and non-protein coding) provided by the ENSEMBL annotations were used. These genes were assigned to phylostratum (ps) 1, and their evolutionary age was assigned a value of 1. On the other end of the spectrum, the number of genes that are common among all four species (ps 4) ranges from 6,902 (pig) to 11,296 (bat). These genes were assigned an evolutionary age of 4. Ps 2 and ps 3 fall between the youngest (ps 1) and oldest (ps 4) genes, and correspond to the branching points of the phylogenetic tree (Fig. 5a). The age of each set of genes was assigned in accordance to its ps. To test if the assignment of homologous genes between any two species was coherent with the expected phylogeny, we performed a hierarchical clustering analysis based on the pair-wise number of homologous genes. We found that the clustering exactly mirrors the phylogeny of the study species (Fig. 5b), supporting the method we used to call homologous genes.
We then calculated the specificity of the transcriptome of each species across development using the transcriptome age index (TAIs) (Fig. 5c). For the forelimbs of mouse and pig, we found that younger genes dominate the transcriptomes of the ridge and paddle stages, while older genes dominate the transcriptome of the bud stage. The forelimb of opossum showed the opposite trend, with older genes dominating the transcriptomes of the ridge and paddle stages and younger genes the transcriptome of the bud stage. Bat demonstrated yet another pattern, with younger genes increasingly dominating the transcriptome as limb development progressed. When we analyzed the statistical significance of the TAIs using VTAI (variance of TAIs) as a test statistic,  we found that the bat, mouse, and pig forelimb trends are highly significant (P-values = 1.11e-4, 1.36e-4, and 3.12e-4, respectively), whereas the opossum forelimb trend is not (P-value = 0.256) (Additional file 6: Figure S5). The same analysis was done for TAIs and VTAI in the hind limbs. In the hind limb, the bat, mouse, and opossum transcriptomes are dominated by older genes at the ridge and paddle stages, and by younger genes at the bud stage. The pig hind limb shows the opposite trend, with the bud stage being dominated by older genes, and the ridge and paddle stages by younger genes. The hind limb trends for bat, mouse, and pig are highly significant (P-values = 5.08e-4, 2.29e-10, and 3.12e-4, respectively) while that of opossum is not significant (P-value = 0.151) (Additional file 7: Figure S6). The lack of statistical significance of opossum results may stem from the reduced number of evolutionary ages that could be assigned to its genes (ps 1 and ps 4). Overall, the genes that are expressed at the paddle stage are younger than those expressed at the ridge stage for the fore- and hind limbs of all species (Fig. 5c).
Similarity of gene expression patterns within (fore- vs. hind limbs) different species – We next determined the degree to which patterns of fore- and hind limb gene expression divergence are similar in the different study species. We found that 4.5% of genes (N = 56) exhibit divergent expression (75th percentile and above) in the fore- and hind limbs of all species during the ridge stage of development, 4.3% of genes (N = 51) exhibit divergent expression during the bud stage, and 8.4% of genes (N = 109) exhibit divergent expression during the paddle stage. These include genes with well-established differences in expression in the fore- and hind limbs (e.g., Tbx4, Tbx5) as well as additional genes with known roles in limb development (e.g., Grem1, Wnt5a) (Additional file 8: Table S1) [15, 61–64]. Of note, these genes do not include Pitx1, a gene known to play a fundamental role in the establishment of fore- and hind limb identity in mammals . The reason for this is that Pitx1 is not present in the pig genome (Sscrofa10.2) we used for our alignment. Therefore, although our results suggest that Pitx1 exhibits divergent expression in the fore- and hind limbs of bat, opossum, and mouse, we have no data on Pitx1 expression in pig. In total, seven genes are divergently expressed in the fore- and hind limbs of all species at all stages (Hand1, Isl1, Myog, Pax1, Tbx4, Tbx5, Tnnt2). This is consistent with literature of some of these genes being differentially expressed in the fore- and hind limb (see Discussion).
Similarity of gene expression patterns among (fore- vs. fore- limb, hind vs. hind limb) different species – We also compared the genes that display the most divergent expression (75th percentile and above) among the limbs of species, and found that several have known roles in limb development. At the forelimb ridge, bud, and paddle stages, 12, 7, and 7 divergently expressed genes have known roles in limb development, respectively (Additional file 9: Table S2). These genes include members of the HoxD family and related genes (Evx2, Lnp), as well as the signaling factor genes Fgf8 and Shh. In the hind limb, the ridge, bud, and paddle stages contain 17, 19, and 18 divergently expressed genes with known roles in limb development, respectively (Additional file 9: Table S2). These genes include additional members of the Hox family and related genes (Hoxa13, Hoxd9, Hoxd11, Hoxd13, Lnp), as well as Fgf8. Two limb genes (Hoxa13 and Med1) are divergently expressed during all three stages of forelimb development, while nine (Ctnnb1, Fbn2, Hoxd11, Med1, Pitx1, Prrx1, Rarb, Tbx4, and Twist1) are divergently expressed at all three stages of hind limb development. Thus, many of the most divergently expressed genes at all stages in the fore- and hindlimbs in our species under study are known to play a role in limb development.
Confirmation of gene expression patterns: To further examine the expression of select genes from the Hoxa and Hoxd clusters (i.e., Hoxd12, Hoxd13, Evx2, Hoxa13), we cloned the coding sequences for bat and opossum and generated species-specific probes for WISH. We focused on the Hoxa and Hoxd clusters because of their well-established role in limb outgrowth and patterning in mice [65, 66]. In general, WISH results were consistent with those of RNA-Seq.
In all species, the Hoxd13 expression domain is confined posteriorly in the fore- and hind limbs before expanding anteriorly. However, the timing and degree of this expansion differs. The degree of anterior-expansion of the Hoxd13 expression domain is greater and more symmetrical in the ridge- and bud-staged forelimbs of opossum (Fig. 6d, e) relative to mouse (Fig. 6a, b), but roughly comparable in paddle-staged opossum (Fig. 6f) and mouse (Fig. 6c) forelimbs. In contrast, the anterior-most boundary of the Hoxd13 expression domain is located more posteriorly in the forelimb paddles of bat (Fig. 6i) relative to mouse (Fig. 6c) and opossum (Fig. 6f). This pattern of early, expanded anterior expression in opossum (not shown) and relatively anteriorly restricted expression in bat are also observed for Hoxa13 (paddle stage shown in Fig. 7). The bat forelimb bud (Fig. 6h) also has a proximodistally wider Hoxd13 expression domain than that of the bat hind limb bud (Fig. 6h’). Hoxd12 (Additional file 10: Figure S10) had similar expression domains in bats, and opossums for stages examined and were generally consistent with published results for mouse Hoxd12 [67–69].
The expression domain of Evx2, a gene that shares regulation with the Hoxd cluster [70–72], is also initially expanded anteriorly in opossum forelimb buds and paddles relative to those of mice (Fig. 8). In mouse, forelimb expression appears in the posterior limb during the early bud stage (Fig. 8a, inset), expands anteriorly during the later bud stage (Fig. 8a), and by the paddle (Fig. 8b) is broadly present across the limb. In contrast, Evx2 is not expressed at ridge (not shown) and early bud stages (Fig. 8a’ inset) of the mouse hind limb, before it appears in the posterior of the limb later in the bud stage (Fig. 8a’) and expands anteriorly in the paddle stage (Fig. 8b’). In opossum, Evx2 is expressed more broadly along the posterior-anterior axis in the fore- (Fig. 8c) than hind limb (Fig. 8c’) bud. Similar to Hoxd13, the opossum Evx2 expression domain is expanded anteriorly and has a more symmetrical domain in the forelimb bud relative to mouse. Furthermore, Evx2 expression is not present in the early hind limb bud stages in mouse (see Fig. 8a’ inset) while it is robustly expressed at this stage in opossum (Fig. 8c’). For additional replicates of Hoxa13, Evx2, and Hoxd13 WISH, see Additional file 11: Figure S7, Additional file 12: Figure S8, Additional file 13: Figure S9.
Similarity of gene expression patterns within (fore- vs. hind limb) and among (fore- vs. fore- limb, hind vs. hind limb) species – We also compared the genes that display the most divergent expression between the fore- and hind limbs of a single species, and among the limbs of species (75th percentile and above). We found that 2.02 to 6.87% of genes exhibit greatly divergent expression in the fore- and hind limb of a single species and among the forelimbs of all species (average = 4.14%). Similarly, 1.79 to 5.35% of genes exhibit greatly divergent expression in the fore- and hind limb of a single species and among the hind limbs of all species (average = 3.29%). At least seven genes with known roles in limb development exhibit greatly divergent expression in the fore- and hind limbs of three of four species and among the limbs of all species (Col2a1, Hoxa13, Mecom, Pitx1, Rarb, Tbx4, Wnt5a, Zbtb16), and at least 5 more in two of four species (Hoxd9, Hoxd13, Msx1, Prrx2, Shh) (Additional file 17: Table S3). As described, we confirmed gene expression patterns of some of these divergently expressed genes (e.g., Hoxa13 and Hoxd13) for some stages and species by in situ hybridization.
In this study we investigated the manner in which expression of the shared genes that regulate limb morphogenesis has been modified during the diversification of mammalian limb morphology. To do this we compared the patterns of gene expression in the fore- and hind limbs of developing bats, mice, opossums, and pigs using RNA-Seq analysis and WISH.
We first investigated the general pattern of divergence of the expression of shared genes among the fore- and hind limbs of a single species. While they vary among species, on the surface our heat-map results suggest that gene expression is generally more conserved at the earliest examined stage of limb development (i.e., ridge) than at later stages. Fore- and hind limb gene expression is most conserved at the ridge stage and diverges at the bud stage in mouse, at the paddle stage in opossum, and appears not to greatly diverge in pig at any examined stage of limb development (i.e., fore- and hind limb of a given developmental stage are more similar to each other than they are to their homologous limbs at other stages). However, while the ridge clusters are statistically significant in mouse, they are not significant in pig or opossum. The lack of significant clustering in opossum at the ridge or any other stage is not that surprising, given that comparable stages of opossum fore- and hind limb development occur at very different stages of overall development when dissimilar organism-wide processes are underway (e.g., skin formation). However, the lack of significant ridge clustering in pig, combined with the significant clustering at the pig bud and paddle stages, suggests that the ridge might not be the most conserved stage of pig development. The pattern in bat, in which the fore- and hind limbs of the paddle stage significantly cluster while those of the ridge and bud do not, also appears to deviate from the early conservation hypothesis. However, it is important to note that forelimb development is often slightly ahead (~12 h in mouse) of hind limb development in placental mammals such as mouse, pig, and bat. Therefore, using a single embryonic stage for the fore- and hind limbs of these species might have somewhat inflated their observed differences in within-species gene expression.
We next investigated the general pattern of divergence of the expression of shared genes among species. Heat-map results suggest that gene expression is conserved among species at the earliest stages of limb development (i.e., ridge), relative to later stages. Specifically, the clustering of species by gene expression during the fore- and hind limb ridge stage more closely follows their phylogenetic relatedness (i.e., bat + pig)  than at later stages. By the paddle stage of limb development, the clustering of species by gene expression more closely follows adult limb morphology, with the highly divergent bat fore- and hind limbs emerging as the outliers. However, while the expression patterns of all species were positively correlated, no specific clusters were significant. This suggests that while the general limb toolkit is conserved among species, the examined species vary in limb gene expression at all stages consistent with their divergent timing of limb development relative to overall development and disparate adult limb morphologies.
Taken together, heat-map results are consistent with gene expression patterns varying more among species during later than earlier stages of limb development, and therefore partially support our first working hypothesis. This part of our first hypothesis is further supported by our calculations of the among-species conservation of gene expression at each developmental stage. These results suggest that for both the fore- and hind limb, the earliest examined stage of limb development (i.e., ridge) displays a highly conserved pattern of gene expression across species that drops significantly as limb development proceeds. However, our heat-map results provide less support for the hypothesis that gene expression patterns vary more within species during later than earlier stages of limb development.
Results of our study of the transcriptomic age index (TAI s ) vary across species and limbs. Of the 8 cases in which we examined transcriptomic age (2 for each species – 1 for the forelimb and 1 for the hind limb), 3 (37.5%; mouse forelimb, pig fore- and hind limb) fit an hourglass model of transcriptomic age, with the middle stage of limb development (i.e., bud) possessing the oldest genes and the youngest (i.e., paddle) and oldest (i.e., ridge) stages the youngest. Four cases (50%; opossum fore- and hind limb, bat hind limb, mouse hind limb) show the opposite pattern, with the middle stage of limb development (i.e., bud) having the youngest genes and the youngest (i.e., paddle) and oldest (i.e., ridge) stages the oldest, and 1 case (12.5%; bat forelimb) displays decreasing transcriptomic age throughout development. However, it is important to note that in every case (100%), transcriptomic age is younger at the latest (i.e., paddle) than earliest (i.e., ridge) stage of limb development.
Results of this study therefore do not support the existence of a phylotypic period of limb development that ends at the onset of chondrogenesis, but do support the hypothesis that gene expression patterns vary more among species during the later than earlier stages of limb development. As a result, these findings are consistent the hypothesis that the hierarchical nature of development translates into increasing variation among species as development progresses [33–39]. However, it is important to note that the temporal range of the proposed phylotypic period remains as controversial as the existence of the period itself ; some authors propose that the phylotypic period may extend throughout the period of organogenesis , while others argue that it may be restricted to the pharyngula stage (~E8.0 in mouse) , or to early somite segmentation (~E8.0 – 8.5 in mouse) , or to the tail-bud stage (~E9.5 – 10.5 in mouse) . While morphological studies of limb development have tended to define the phylotypic period more broadly [43, 44], recent embryo-wide molecular studies in vertebrates suggest that the phylotypic period might actually be temporally restricted to the earliest stages of limb development (E8.0-E8.5) [51, 73]. At least the among species results of this study would be consistent with the existence of a phylotypic period for the limb that ends by onset of the bud stage, or roughly E10.5 in mouse.
Going beyond general patterns of gene expression, we also investigated the relative expression levels of specific genes in the fore- and hind limbs of single species and among the limbs of different species. Although we identified some limb-patterning genes (e.g., Hand1, Pax1, Tbx4, Tbx5) that are differentially expressed in the fore- and hind limbs of all examined species, these genes represent a small subset (0.5 to 8.4%) of the differentially expressed genes in all species. Some of these identified genes were well-supported in the literature. For example, a DGE-tag study of Myotis ricketti (Rickett’s big-footed bat) found Tbx5 to be expressed more highly in the forelimb digits of fetal bats, while Tbx4 was more highly expressed in the hindlimb digits . Tbx5 is well established to be forelimb restricted, while Tbx4 is hindlimb restricted in mice, opossums, and chickens [23, 63, 64], providing support for our RNA-Seq results. We also identified a suite of key limb-patterning genes that are differentially expressed among the homologous limbs of all species. Therefore, results of this study are consistent with the hypothesis that the expression patterns of several limb-patterning genes significantly differ between the fore- and hind limbs of a single mammalian species and among the homologous limbs of different species.
The results of the RNA-Seq assays were generally supported by WISH from this study and previous studies. Results of this study confirmed that several Hoxa and Hoxd cluster genes exhibit distinct expression domains among species. Previous studies demonstrated that Prrx1 (also known as MHox or Prx1) is expressed in similar domains during early stages of bat and mouse limb development, but that its’ expression expands in the hand-plate of bats relative to those of mice . Fgf8 is normally expressed in the AER (Apical Ectodermal Ridge) of mammalian limbs [78, 79]. Relative to mouse, Fgf8 expression in pig has been shown to be reduced at later stages of limb development . In contrast, bat forelimbs have previously been shown to have a wider Fgf8 expression domain in the AER relative to the forelimbs of mice at similar stages .
Results also suggest that only a small subset of genes (1.79 to 6.87%) displays divergent expression within (e.g., between the fore- and hind limb of a single species) and among (e.g., between the forelimbs of two different species) species. This finding suggests that the genes that display divergent expression in the fore- and hind limb of a single species are not likely to also display divergent expression across the homologous limbs of multiple species. Therefore, in contrast to our working hypothesis, results of this study are consistent with a scenario in which evolutionary divergence of limb form within (i.e., of the fore- and hind limb) and among species (i.e., of the forelimbs) is driven by changes in the expression of different genes. Taken together, results of this study are consistent with a scenario in which a small subset of limb genes controls fore- vs. hind limb identity, and a small but different subset controls the evolution of limb morphology across multiple species. This result, especially when combined with the heat-map results described above, suggests that different processes might be controlling variation in limb gene expression within- and among- species.
While these results are intriguing, one potential caveat of this study is that we were forced to use heterogeneous methods to align the RNA-seq reads of our study species. We first attempted to align reads of all species directly to reference genomes. While this worked well (i.e., high alignment rate) for mouse, pig, and opossum, species for which high quality genomes are available, this did not work well for bat. As the bat species we used, Carollia perspicillata, does not have a published reference genome, we initially tried to align to the published reference genome for Myotis lucifugus, a phylogenetically-related bat species. However, only 7-9% of Carollia reads aligned with the Myotis genome (Additional file 14: Table S7). To overcome this, we generated a de novo transcriptome for Carollia, which increased the alignment rate to a much more reasonable ~80% for the bat reads, and used these data for subsequent analyses. While the use of different mapping techniques for bat versus pig, mouse, and opossum could potentially introduce bias into our analysis, previous studies suggest that the aligning methods we used under the conditions we used them do not result in noticeable differences . Furthermore, if the difference in methods was introducing significant bias, then we would expect for bat reads to fall as outliers in all cross-species comparisons. Instead, we found that bat reads often clustered with pig or mouse to the exclusion of other species. Therefore, while it is important to note that we were forced to use heterogeneous methods to analyze our data, our data are not consistent with these different methods biasing our results. For further discussion, see Methods and Additional file 1.
Using RNA-Seq, we generated transcriptomes for the ridge, bud, and paddle stages of limb development in mouse, opossum, bat, and pig to explore the hypothesis that mammalian limb diversification has occurred through the differential expression of shared conserved genes. Generally, there is greater variation among species at later stages (paddle) of development than at earlier stages (ridge). In addition, genes expressed at later stages tend to be younger in evolutionary age than those expressed at earlier (ridge) stages. Several key limb-patterning genes are differentially expressed among the homologous limbs of the species under study, though only a small subset of these shared genes are differentially expressed in the fore- and hind limbs of all species. In addition, only a small subset of these shared genes are differentially expressed in the fore- and hind limbs of the same species. Our results generally support the hypothesis that variation among species increases as development progresses, but do not support the presence of a phylotypic period of limb development that ends at chondrogenesis.
Embryo collection and staging – All procedures were within the guidelines of the University of Illinois Institutional Animal Care and Use Committee (IACUC). Embryonic mice (ICR strain, Taconic) and opossums were obtained from timed matings [81, 82] in breeding colonies housed in the Sears Lab at the University of Illinois. Pig embryos were obtained through timed inseminations following ovulations at the University of Illinois pig farm . For pig, mouse, and opossum, embryos were dissected out and placed in RNALater at 4 °C overnight before being frozen until the limbs were dissected off. Bat embryos were obtained from field collections in Trinidad and staged according to Cretekos et al. 2005 . Forelimbs and hindlimbs were dissected and placed in RNALater (Qiagen) and stored at room temperature. Limbs from the ridge (Wanek stage 2), bud (Wanek stage 3/4), and paddle (Wanek stage 6) stages of development were obtained for all species . In bat these stages fall around Carollia stages (CS) 13, 14, and 15 for both the fore- and hind limb, in mouse around embryonic days (E) 10, 10.5, and 11.5 for both the fore- and hind limb, in opossum around stages 27, 28, and 29 for the forelimb, and 30, 31, and 32 for the hind limb, and in pig around embryonic days (E) 21.5, 23.5, and 25.5 for both the fore- and hind limb [57, 84–87]. At least three biological replicates were collected for each species/stage, though some were excluded from analysis post-library construction. Forelimbs from a single embryo and hindlimbs from a single embryo were combined for the bud and paddle stages. For ridge stages, limbs from 2 embryos were pooled.
Embryos for whole-mount in situ hybridization (WISH) were collected at the appropriate stages, fixed in 4% paraformaldehyde (PFA) overnight and dehydrated through a methanol series before being processed for WISH (see below).
RNA sequencing - Limbs were removed from embryos as above and stored in RNALater (Life Technologies) at −20 °C until further processing. RNA was extracted from tissues using the E.Z.N.A. Total RNA Kit I (OMEGA bio-tek #R6834), and converted into RNASeq libraries with the Illumina TruSeq RNA Sample Preparation Kit (Illumina RS-122-2001). Before RNASeq library construction, a few RNA samples were checked for quality (Bioanalyzer). Libraries were sequenced on an Illumina HiSeq 2500 housed in the Roy G. Carver Biotechnology Center at the University of Illinois under the supervision of Dr. Alvaro Hernandez. For all species, the resulting single-end RNASeq sequences were pre-processed to remove Illumina adaptors and low quality bases (score <20) from the 3’ read end. The pre-processed libraries where then aligned to their corresponding reference genomes. For mouse, opossum, and pig, we used the Ensembl reference genomes and annotation files corresponding to assemblies, GRCm38, BROADO5, and Sscrofa10.2, respectively. As the bat species we are using (Carollia perspicillata) does not have a published genome, we initially used Myotis lucifugus (The Little Brown Bat, whose genome was sequenced ) as a reference genome (Myoluc2.0), but the alignment rate with TopHat  for Carollia sequences was very low (~7-9% for bud and paddle stages). Therefore, we next performed a de novo transcriptome assembly with Carollia reads. To do this, we pooled all bat libraries and fed them into Trinity , a de-novo assembly tool. Trinity has been used by other groups to assemble bat RNASeq transcriptomes from various tissues , and several bat mitochondrial genomes (. Trinity has become increasingly popular for other species for which a genome is not available . Pooled cleaned libraries were combined into a single file, and duplicated sequences removed to reduce computational requirements. Trinity predicted 350,733 gene transcripts which were filtered to keep only those matching the protein sequences of the SwissProt database (blastx, E-value < 1e-20) . This resulted in 88,930 gene isoforms. These were used as a reference to map the RNA-Seq libraries. For bat, read alignment and gene expression were computed with RSEM  using as reference sequences the 88,930 gene isoforms previously described. RNA-Seq data has been deposited in the NCBI Gene Expression Omnibus (GEO) under accession number GSE71390 .
For species with a reference genome (mouse, opossum, and pig) we aligned the reads using STAR (Spliced Transcripts Alignment to a Reference)  and then used Cufflinks  to compute their gene expression. For all four species, gene expression values were normalized as fragments per kilobase of transcript per million mapped reads (FPKM). To reduce noise when quantifying FPKM values we used all reads aligned to genes, and did not discard reads aligned to non-constitutive exons or non-orthologous gene fragments. To account for differences in the mass composition between samples we used DESeq normalization as described by . All analyses were performed in R .
The aligner used for mouse, pig, and opossum (STAR), could not be used for the bat species Carollia because a species-specific reference genome for Carollia is not available. Instead, BOWTIE, which is internally used by RSEM, was used. A comparative analysis of various aligners, including STAR and TOPHAT (which, like RSEM, relies on BOWTIE for alignment), shows that noticeable differences between aligners only result when reads are mapped to unannotated junctions . In our methods we did not use unannotated junctions (bat reads were aligned against the de novo transcriptome, and all other species were aligned only to annotated reference genes). Thus, using BOWTIE for bat and STAR for all other species does not introduce a significant source of bias in our analyses.
Analyses, Divergence of the expression pattern of shared genes within and among species - To identify when the expression pattern of shared genes diverges among the fore- and hind limbs of a single species and among the limbs of different species, we first graphically visualized the similarity in expression of a set of 6,583 genes (orthologous to all four species) in all four species (bat, mouse, opossum, pig) at each stage of limb development (ridge, bud, paddle) using a series of heat-map analyses, where distances between any to gene expression profiles were computed using Spearman correlation coefficients. Testing for the significance of the clustering was done using pvclust . To determine these orthologous genes, the transcriptome from Carollia was used as a reference to align the transcriptomes of mouse, pig, and opossum (blastn, E-value 1e-20). Bat genes matching more than one orthologous gene in any other species were filtered out, and then the bat genes that have orthologous sequences in all the other species determined. This gave a set of 6,583 genes found in all 4 species.
We then used the set of 6,583 genes (orthologous to all four species) to calculate the among-species conservation of gene expression at each developmental stage, using the mean of all species pairwise Spearman coefficients (c):
Where r i, j is the Spearman coefficient between species i and j at a given stage, and
k is the total number of species under study. We used the Spearman rather than Pearson coefficient, as the former is more robust against outliers. To test the robustness of any differences in Spearman coefficients among species with respect to the selection of orthologous genes, we randomly sub-sampled 500 sets of genes at early and late stages at intensities ranging from 50 to 100% of all genes.
To investigate the divergence of the expression pattern of shared genes among species, we computed the evolutionary age of the genes of each species. To do this we determined the most distant phylogenetic node among bat, pig, mouse, and opossum that contained a detectable homologue for each gene. FASTA sequences of all species genes were pooled to form a single database to which each species genome was aligned (blast hit, E-value 1e-5 for homologue detection) . In this analysis, youth corresponds to genes that have more recently evolved and, therefore, are more likely to be specific to a single species. Conversely, old genes are those that have been inherited from a common ancestor and are present among all descendant species. The specificity of the transcriptome of each species was then quantified across limb development (ridge, bud, paddle) using the transcriptome age index (TAI), which is the sum of each gene evolutionary age weighted by its expression . As a result, a smaller TAI value corresponds to a younger age, and a larger TAI value an older age. For a given species at a given stage, s, TAI s is mathematically defined as:
Here, ps i and e i are the age and expression values of gene i, respectively. The total number of genes is represented by n. We then analyzed the statistical significant of the TAIs trend of each species using the procedure proposed by , in which the variance of TAIs across stages (ridge, budge, and paddle), VTAI, is used as a test statistic. The null distribution for this analysis was obtained by sampling 1000 VTAI surrogates. Each surrogate was generated by randomly permuting the ps assignations. The null distribution was modeled as a gamma distribution, and its parameters estimated using the MASS library in R .
Similarity of gene expression patterns within species – To determine the degree to which patterns of fore- and hind limb gene expression divergence are similar between species we identified those genes with the most divergent expression (75% percentile and above) between the fore- and hind limb for each species. We then performed a series of pairwise comparisons in which we compared the resulting lists of genes for two species and determined the percentage of the total number of genes that are present in both.
Similarity of gene expression patterns within and among species – We also determined the degree to which patterns of gene expression divergence are similar within (e.g., between the fore- and hind limb of a single species) and among (e.g., between the forelimbs of two different species) species. To do this we identified those genes with the most divergent expression (75% percentile and above): between the fore- and hind limb for each species (within species, as described above), among the forelimbs of all species (among species), and among the hind limbs of all species (among species). We then compared the within and among species lists and determined the percentage of the total number of genes that are present in both. We also used the DAVID Ontology database (http://david.abcc.ncifcrf.gov/) to identify a subset of genes from the lists with known roles in limb development [99, 100].
Whole-mount in situ hybridization (WISH)
To confirm the RNA-Seq expression data of some of our divergently expressed candidate genes, WISH was performed on embryos at the ridge, bud and paddle stages for the species investigated (See Additional file 15: Table S6). Often, mouse differs enough from bat and opossum that it is necessary to make species-specific probes for WISH. Coding sequences were found in the NCBI Nucleotide database and primers designed specific to each species using NCBI PrimerBLAST . (See Additional file 16: Table S4 for accession numbers and primer sequences). The primers were synthesized by Sigma-Aldrich. RNA was purified from bat and opossum limbs and used to generate cDNA using the SuperScript III Reverse Transcriptase kit (Invitrogen) following the manufacturer’s instructions. PCR was done to amplify the cDNA for the gene of interest using the species-specific primers and standard PCR methods. Genes were cloned into pGem T easy (see Additional file 16: Table S4 for primer sequences and accession numbers) and sequenced. After sequencing the direction and sequence identity were confirmed using NCBI Blast. Next, plasmid DNA was linearized with the appropriate restriction enzyme (NotI or SpeI) and an in vitro transcription reaction to generate antisense mRNA probes was performed using Roche or Promega reagents. Probes were all digoxigenin labeled. After synthesis, probes were purified by ethanol precipitation and checked on a NanoDrop for concentration and RNA Quality. Probes were stored at −80 °C until use.
To perform the in situ hybridization, standard methods were based on the following protocols [102–105]. BM Purple (Roche) was used to develop the reactions. After development (assessed by purple/blue staining where the probe has bound), the embryos were photographed on a Leica camera microscope and fixed in 4% PFA for long-term storage.
Apical ectodermal ridge
Database for annotation, visualization, and integrated discovery
Digital gene expression tag
Fragments per Kilobase of transcript per million mapped reads
Modern applied statistics with S
National center for biotechnology information
Polymerase chain reaction
RNA-Seq by Expectation Maximization
Spliced transcripts alignment to a reference
Transcriptome age index
Variance of TAIs
Whole mount in situ hybridization
Shubin N, Tabin C, Carroll S. Deep homology and the origins of evolutionary novelty. Nature. 2009;457(7231):818–23.
Hockman D, Cretekos CJ, Mason MK, Behringer RR, Jacobs DS, Illing N. A second wave of sonic hedgehog expression during the development of the bat limb. Proc Natl Acad Sci U S A. 2008;105(44):16982–7.
Sears KE, Behringer RR, Rasweiler JJ, Niswander LA. Development of bat flight: morphologic and molecular evolution of bat wing digits. Proc Natl Acad Sci U S A. 2006;103(17):6581–6.
Thewissen JGM, Cohn MJ, Stevens ME, Bajpai S, Heyning J, Horton WEJ. Developmental basis for hind-limb loss in dolphins and origin of the cetacean bodyplan. Proc Natl Acad Sci U S A. 2006;103(22):8414–8.
Carroll SB. Evolution at two levels: on genes and form. PLoS Biol. 2005;3(7):1159–66.
King MC, Wilson AC. Evolution at two levels in humans and chimpanzees. Science. 1975;188(4184):107–16.
Gompel N, Prud’homme B, Wittkopp PJ, Kassner VA, Carroll SB. Chance caught on the wing: cis-regulatory evolution and the origin of pigment patterns in Drosophila. Nature. 2005;433:481–7.
Cooper KL, Sears KE, Uygur A, Maier J, Baczkowski KS, Brosnahan M, Antczak D, Skidmore JA, Tabin CJ. Patterning and post-patterning modes of evolutionary digit loss in mammals. Nature. 2014;511(7507):41–5.
Tschopp P, Sherratt E, Sanger TJ, Groner AC, Aspiras AC, Hu JK, Pourquie O, Gros J, Tabin CJ. A relative shift in cloacal location repositions external genitalia in amniote evolution. Nature. 2014;516(7531):391–4.
Gehrke AR, Schneider I, de la Calle-Mustienes E, Tena JJ, Gomez-Marin C, Chandran M, Nakamura T, Braasch I, Postlethwait JH, Gomez-Skarmeta JL, et al. Deep conservation of wrist and digit enhancers in fish. Proc Natl Acad Sci U S A. 2015;112(3):803–8.
Cretekos CJ, Wang Y, Green ED, Martin JF, Rasweiler JJ, Behringer RR. Regulatory divergence modifies limb length between mammals. Genes Dev. 2008;22(2):141–51.
Cretekos CJ, Deng JM, Green ED, Program NCS, Rasweiler JJ, Behringer RR. Isolation, genomic structure and developmental expression of Fgf8 in the short-tailed fruit bat, Carollia perspicillata. Int J Dev Biol. 2007;51(4):333–8.
Polly PD. Limbs in mammalian evolution. In: Hall BK, editor. Fins into limbs: evolution, development, and transformation. Chicago: University of Chicago Press; 2007. p. 245–68.
Zeller R. The temporal dynamics of vertebrate limb development, teratogenesis and evolution. Curr Opin Genet Dev. 2010;20(4):384–90.
Zeller R, Lopez-Rios J, Zuniga A. Vertebrate limb bud development: moving towards integrative analysis of organogenesis. Nat Rev Genet. 2009;10(12):845–58.
Ray R, Capecchi MR. An examination of the chiropteran HoxD locus from an evolutionary perspective. Evol Dev. 2008;10(6):657–70.
Weatherbee SD, Behringer RR, Rasweiler JJ, Niswander L. Interdigital webbing retention in bat wings illustrates genetic changes underlying amniote limb diversification. Proc Natl Acad Sci U S A. 2006;103(41):15103–7.
Sears KE. Molecular determinants of bat wing development. Cells Tissues Organs. 2008;187:6–12.
Chen CH, Cretekos CJ, Rasweiler JJ, Behringer RR. Hoxd13 expression in the developing limbs of the short-tailed fruit bat, Carollia perspicillata. Evol Dev. 2005;7(2):130–41.
Doroba CK, Sears KE. The divergent developmental of the apical ectodermal ridge in the marsupial Monodelphis domestica. Anat Rec. 2010;293(8):1325–32.
Hubler M, Molineaux AC, Keyte A, Schecker T, Sears KE. Development of the marsupial shoulder girdle complex: a case study in Monodelphis domestica. Evol Dev. 2013;15(1):18–27.
Sears KE, Patel A, Hübler M, Cao X, VandeBerg JL, Zhong S. Disparate Igf1 expression and growth in the fore- and hind limbs of a marsupial (Monodelphis domestica). J Exp Zool B: Mol Dev Evol. 2012;318(4):279–93.
Keyte AL, Smith KK. Developmental origins of precocial forelimbs in marsupial neonates. Development. 2010;137(24):4283–94.
Wang Z, Dai M, Wang Y, Cooper KL, Zhu T, Dong D, Zhang J, Zhang S. Unique expression patterns of multiple key genes associated with the evolution of mammalian flight. Proc Biol Sci/R Soc. 2014;281(1783):20133133.
Taher L, Collette NM, Murugesh D, Maxwell E, Ovcharenko I, Loots GG. Global gene expression analysis of murine limb development. PLoS One. 2011;6(12):e28358.
Sears KE, Doroba CK, Xie D, Zhong S. Molecular determinants of marsupial limb integration and constraint. In: Müller J, Asher R, editors. From clone to bone: The synergy of morphological and molecular tools in paleobiology. Cambridge: Cambridge University Press; 2012. p. 257–78.
Eckalbar WL, Schlebusch SA, Mason MK, Gill Z, Parker AV, Booker BM, Nishizaki S, Muswamba-Nday C, Terhune E, Nevonen KA, et al. Transcriptomic and epigenomic characterization of the developing bat wing. Nat Genet. 2016;48(5):528–36.
Jenkins FA, Parrington FR. The postcranial skeletons of the Triassic mammals eozostrodon, megazostrodon, and erythrotherium. Philos Trans R Soc Lond B. 1976;273(926):387–431.
Kemp TS. The origin and evolution of mammals. New York: Oxford University Press; 2005.
Lopez-Rios J, Duchesne A, Speziale D, Andrey G, Peterson KA, Germann P, Unal E, Liu J, Floriot S, Barbey S, et al. Attenuated sensing of SHH by Ptch1 underlies evolution of bovine limbs. Nature. 2014;511(7507):46–51.
Cooper LN, Cretekos CC, Sears KE. The evo-devo of mammalian flight. WIREs Dev Biol. 2012;1(5):773–9.
Stopper GF, Wagner GP. Of chicken wings and frog legs: a smorgasbord of evolutionary variation in mechanisms of tetrapod limb development. Dev Biol. 2005;288:21–39.
von Baer KE. Entwicklungsgeschichte der Thiere: Beobachtung und Reflexion. Konigsberg: Borntrager; 1828.
Garstang W. The theory of recapitulation: a critical restatement of the Biogenetic law. J Exp Zool. 1922;291:195–204.
de Beer GR. Embryos and ancestors, revised edition. Oxford: Oxford University Press; 1954.
Wolpert L. Positional information and pattern formation. Philos Trans R Soc Lond Ser B Biol Sci. 1981;295(1078):441–50.
Raff RA. The shape of life: genes, development, and the evolution of animal form. Chicago: University of Chicago Press; 1996.
Kalinka AT, Tomancak P. The evolution of early animal embryos: conservation or divergence? Trends Ecol Evol. 2012;27(7):385–93.
Poe S. Test of Von Baer’s law of the conservation of early development. Evolution. 2006;60(11):2239–45.
Roux J, Robinson-Rechavi M. Developmental constraints on vertebrate genome evolution. PLoS Genet. 2008;4(12):e1000311.
Richardson MK. A phylotypic stage for all animals? Dev Cell. 2012;22(5):903–4.
Duboule D. Temporal colinearity and the phylotypic progression: a basis for the stability of a vertebrate Bauplan and the evolution of morphologies through heterochrony. Development. 1994;1994:135–42.
Galis F, Jacques JMA, Metz JAJ. Why five fingers? Evolutionary constraints on digit numbers. Trends Ecol Evol. 2001;16(11):637–46.
Galis F, Metz JAJ. Testing the vulnerability of the phylotypic stage: on modularity and evolutionary conservation. J Exp Zool (Mol Dev Evol). 2001;291:195–204.
Hazkani-Covo E, Wool D, Graur D. In search of the vertebrate phylotypic stage: a molecular examination of the developmental hourglass model and von Baer’s third law. J Exp Zool B. 2005;304:150–8.
Galis F, Alphen JJM, Metz JAJ. Digit reduction: via repatterning or developmental arrest. Evol Dev. 2002;4(4):249–51.
Quint M, Drost HG, Gabel A, Ullrich KK, Bonn M, Grosse I. A transcriptomic hourglass in plant embryogenesis. Nature. 2012;490(7418):98–101.
Kalinka AT, Varga KM, Gerrard DT, Preibisch S, Corcoran DL, Jarrells J, Ohler U, Bergman CM, Tomancak P. Gene expression divergence recapitulates the developmental hourglass model. Nature. 2010;468(7325):811–4.
Domazet-Loso T, Tautz D. A phylogenetically based transcriptome age index mirrors ontogenetic divergence patterns. Nature. 2010;468(7325):815–8.
Levin M, Hashimshony T, Wagner F, Yanai I. Developmental milestones punctuate gene expression in the Caenorhabditis embryo. Dev Cell. 2012;22(5):1101–8.
Irie N, Kuratani S. Comparative transcriptome analysis reveals vertebrate phylotypic period during organogenesis. Nat Commun. 2011;2:248.
Sears KE. Differences in the timing of prechondrogenic limb development in mammals: the marsupial-placental dichotomy resolved. Evolution. 2009;63(8):2193–200.
Sears KE, Maier JA, Rivas-Astroza M, Poe R, Zhong S, Kosog K, Marcot JD, Behringer RR, Cretekos CJ, Rasweiler JJ, et al. The relationship between gene network structure and expression variation among individuals and species. PLoS Genet. 2015;11(8):e1005398.
Ross D, Marcot JD, Betteridge KJ, Nascone-Yoder N, Bailey CS, Sears KE. Constraints on mammalian forelimb development: Insights from developmental disparity. Evolution. 2013;67(12):3645–52.
Sears KE. Novel insights into the regulation of limb development from “natural” mammalian mutants. BioEssays. 2011;33(5):327–31.
Wanek N, Muneoka K, Holler-Dinsmore G, Burton R, Bryant SV. A staging system for mouse limb development. J Exp Zool. 1989;249(1):41–9.
Cretekos CJ, Weatherbee SD, Chen CH, Badwaik NK, Niswander L, Behringer RR, Rasweiler JJ. Embryonic staging system for the short-tailed fruit bat, Carollia perspicillata, a model organism for the mammalian order Chiroptera, based upon timed pregnancies in captive-bred animals. Dev Dyn. 2005;233(3):721–38.
Sears KE, Bormet AK, Rockwell A, Powers LE, Noelle Cooper L, Wheeler MB. Developmental basis of mammalian digit reduction: a case study in pigs. Evol Dev. 2011;13(6):533–41.
Suzuki R, Shimodaira H. Pvclust: an R package for assessing the uncertainty in hierarchical clustering. Bioinformatics. 2006;22(12):1540–2.
Meredith RW, Janecka JE, Gatesy J, Ryder OA, Fisher CA, Teeling EC, Goodbla A, Eizirik E, Simao TL, Stadler T, et al. Impacts of the Cretaceous Terrestrial Revolution and KPg extinction on mammal diversification. Science. 2011;334(6055):521–4.
Verheyden JM, Sun X. An Fgf/Gremlin inhibitory feedback loop triggers termination of limb bud outgrowth. Nature. 2008;454(7204):638–41.
Gros J, Hu JK, Vinegoni C, Feruglio PF, Weissleder R, Tabin CJ. WNT5A/JNK and FGF/MAPK pathways regulate the cellular events shaping the vertebrate limb bud. Curr Biol. 2010;20(22):1993–2002.
Minguillon C, Del Buono J, Logan MP. Tbx5 and Tbx4 are not sufficient to determine limb-specific morphologies but have common roles in initiating limb outgrowth. Dev Cell. 2005;8(1):75–84.
Rodriguez-Esteban C, Tsukui T, Yonei S, Magallon J, Tamura K, Izpisua Belmonte JC. The T-box genes Tbx4 and Tbx5 regulate limb outgrowth and identity. Nature. 1999;398(6730):814–8.
Graham A, Papalopulu N, Krumlauf R. The murine and Drosophila homeobox gene complexes have common features of organization and expression. Cell. 1989;57(3):367–78.
Montavon T, Soshnikova N. Hox gene regulation and timing in embryogenesis. Semin Cell Dev Biol. 2014;34:76–84.
Albrecht AN, Schwabe GC, Stricker S, Boddrich A, Wanker EE, Mundlos S. The synpolydactyly homolog (spdh) mutation in the mouse -- a defect in patterning and growth of limb cartilage elements. Mech Dev. 2002;112(1–2):53–67.
Chen Y, Knezevic V, Ervin V, Hutson R, Ward Y, Mackem S. Direct interaction with Hoxd proteins reverses Gli3-repressor function to promote digit formation downstream of Shh. Development. 2004;131(10):2339–47.
Knezevic V, De Santo R, Schughart K, Huffstadt U, Chiang C, Mahon KA, Mackem S. Hoxd-12 differentially affects preaxial and postaxial chondrogenic branches in the limb and regulates Sonic hedgehog in a positive feedback loop. Development. 1997;124(22):4523–36.
Kmita M, Fraudeau N, Herault Y, Duboule D. Serial deletions and duplications suggest a mechanism for the collinearity of Hoxd genes in limbs. Nature. 2002;420(6912):145–50.
Lehoczky JA, Williams ME, Innis JW. Conserved expression domains for genes upstream and within the HoxA and HoxD clusters suggests a long-range enhancer existed before cluster duplication. Evol Dev. 2004;6(6):423–30.
Spitz F, Duboule D. Global control regions and regulatory landscapes in vertebrate development and evolution. Adv Genet. 2008;61:175–205.
Irie N, Sehara-Fujisawa A. The vertebrate phylotypic stage and an early bilaterian-related stage in mouse embryogenesis defined by genomic information. BMC Biol. 2007;5(1):1–8.
Ballard WB. Morphogenetic movements and fate maps of vertebrates. Am Zool. 1981;21:391–9.
Wolpert L. The triumph of the embryo. Oxford England. New York: Oxford University Press; 1991.
Slack JM, Holland PW, Graham CF. The zootype and the phylotypic stage. Nature. 1993;361(6412):490–2.
Wang Z, Dong D, Ru B, Young RL, Han N, Guo T, Zhang S. Digital gene expression tag profiling of bat digits provides robust candidates contributing to wing formation. BMC Genomics. 2010;11:619.
Pajni-Underwood S, Wilson CP, Elder C, Mishina Y, Lewandoski M. BMP signals control limb bud interdigital programmed cell death by regulating FGF signaling. Development. 2007;134(12):2359–68.
Molineaux AC, Maier JA, Schecker T, Sears KE. Exogenous retinoic acid induces digit reduction in opossums (Monodelphis domestica) by disrupting cell death and proliferation, and apical ectodermal ridge and zone of polarizing activity function. Birth Defects Res A, Clin Mol Teratol. 2015;103(3):225–34.
Haas BJ, Papanicolaou A, Yassour M, Grabherr M, Blood PD, Bowden J, Couger MB, Eccles D, Li B, Lieber M, et al. De novo transcript sequence reconstruction from RNA-seq using the Trinity platform for reference generation and analysis. Nat Protoc. 2013;8(8):1494–512.
Keyte AL, Smith KK. Basic maintenance and breeding of the opossum monodelphis domestica. CSH Protoc. 2008;2008:pdb prot5073.
Suckow MA, Danneman P, Brayton C. The laboratory mouse. Boca Raton: CRC Press; 2001.
Artificial Insemination (AI) [http://livestocktrail.illinois.edu/swinerepronet/paperDisplay.cfm?ContentID=6267]. Accessed May-June 2016.
Patten BM. Embryology of the Pig. CoPhiladelphia, USA: Blakiston; 1931.
McCrady E. The embryology of the Opossum. Philadelphia: Wistar Institute of Anatomy and Biology; 1938.
Mate KE, Robinson ES, Vandeberg JL, Pedersen RA. Timetable of in vivo embryonic development in the grey short-tailed opossum (Monodelphis domestica). Mol Reprod Dev. 1994;39(4):365–74.
Theiler K. The house mouse. New York: Springer; 1972.
Lindblad-Toh K, Garber M, Zuk O, Lin MF, Parker BJ, Washietl S, Kheradpour P, Ernst J, Jordan G, Mauceli E, et al. A high-resolution map of human evolutionary constraint using 29 mammals. Nature. 2011;478(7370):476–82.
Trapnell C, Roberts A, Goff L, Pertea G, Kim D, Kelley DR, Pimentel H, Salzberg SL, Rinn JL, Pachter L. Differential gene and transcript expression analysis of RNA-seq experiments with TopHat and Cufflinks. Nat Protoc. 2012;7(3):562–78.
Grabherr MG, Haas BJ, Yassour M, Levin JZ, Thompson DA, Amit I, Adiconis X, Fan L, Raychowdhury R, Zeng Q, et al. Full-length transcriptome assembly from RNA-Seq data without a reference genome. Nat Biotechnol. 2011;29(7):644–52.
Lee AK, Kulcsar KA, Elliott O, Khiabanian H, Nagle ER, Jones ME, Amman BR, Sanchez-Lockhart M, Towner JS, Palacios G, et al. De novo transcriptome reconstruction and annotation of the Egyptian rousette bat. BMC Genomics. 2015;16:1033.
Shi H, Zhang S, Mao X: The complete mitochondrial genome of the king horseshoe bat (Rhinolophus rex) using next-generation sequencing and Sanger sequencing. Mitochondrial DNA 2015:1–2.
Bairoch A, Apweiler R. The SWISS-PROT protein sequence data bank and its supplement TrEMBL. Nucleic Acids Res. 1997;25(1):31–6.
Li B, Dewey CN. RSEM: accurate transcript quantification from RNA-Seq data with or without a reference genome. BMC Bioinformatics. 2011;12:323.
Dobin A, Davis CA, Schlesinger F, Drenkow J, Zaleski C, Jha S, Batut P, Chaisson M, Gingeras TR. STAR: ultrafast universal RNA-seq aligner. Bioinformatics. 2013;29(1):15–21.
Anders S, Huber W. Differential expression analysis for sequence count data. Genome Biol. 2010;11(10):R106.
Team RDC. R: a language and environment for statistical computing. Vienna, Austria: R Foundation for Statistical Computing; 2011.
Venables WN, Ripley BD, Venables WN, SpringerLink (Online service). Modern applied statistics with S. In: Statistics and computing. 4th ed. New York: Springer; 2002. 1 online resource (xi, 495 pages).
Huang DW, Sherman BT, Lempicki RA. Systematic and integrative analysis of large gene lists using DAVID bioinformatics resources. Nat Protoc. 2009;4(1):44–57.
Huang DW, Sherman BT, Lempicki RA. Bioinformatics enrichment tools: paths toward the comprehensive functional analysis of large gene lists. Nucleic Acids Res. 2009;37(1):1–13.
Ye J, Coulouris G, Zaretskaya I, Cutcutache I, Rozen S, Madden TL. Primer-BLAST: a tool to design target-specific primers for polymerase chain reaction. BMC Bioinformatics. 2012;13:134.
Keyte AL, Smith KK. Whole-mount in situ hybridization in monodelphis embryos. CSH Protoc. 2008;2008:pdb prot5076.
Murtaugh LC, Chyung JH, Lassar AB. Sonic hedgehog promotes somitic chondrogenesis by altering the cellular response to BMP signaling. Genes Dev. 1999;13(2):225–37.
Nieto MA, Patel K, Wilkinson DG. In situ hybridization analysis of chick embryos in whole mount and tissue sections. Methods Cell Biol. 1996;51:219–35.
Rasweiler JJ, Cretekos CJ, Behringer RR. Whole-mount in situ hybridization of short-tailed fruit bat (Carollia perspicillata) embryos with RNA probes. Cold Spring Harb Protoc. 2009;2009(3):pdb prot5164.
We thank Simeon Williams for field support during bat collections and the Wildlife Section, Forestry Division, Ministry of Agriculture, Land, and Marine Resources (currently in the Ministry of Public Utilities and the Environment) of the Republic of Trinidad and Tobago for the issuance of required bat collecting and export permits. We also thank Lori Rund and Jonathan Mosley for assistance with pig breeding. For WISH probes, we thank Dr. Brian Harfe for mouse Hoxd13 and Dr. Denis Duboule for mouse Evx2 plasmids.
This work was supported by the Division of Integrative Organismal Systems of the National Science Foundation (grant number 1257873 to K.E.S. and S.Z., 1256423 to R.R.B., and 1255926 to C.J.C.).
Availability of data and materials
RNA-Seq data from these experiments has been deposited in the GEO database under accession number GSE71390. Additional data and code will be deposited in DRYAD/NCBI databases and additional libraries to GEO. Requests for plasmids can be sent to KES. Additional files can be found in Additional file 1, Additional file 2: Figures S1, Additional file 3: Figure S2, Additional file 4: Figure S3, Additional file 5: Figure S4, Additional file 6: Figure S5, Additional file 7: Figure S6, and Additional file 8: Tables S1, Additional file 9: Tables S2, Additional file 17: Tables S3, Additional file 16: Tables S4, Additional file 18: Tables S5, Additional file 15: Tables S6.
RRB, CJC, and JJR collected bat embryos; JAM and AD bred and collected mouse and opossum embryos; JAM and KES purified RNA and constructed RNA-Seq libraries; MR and XC performed bioinformatics analysis of RNA-Seq data; JAM, JD, AD, and PO cloned coding sequences and generated probes; JAM and PO performed WISH; KES, RRB, SZ, CJC, and RJR conceived and designed the experiments, KES, JAM, MR, XC, RRB, and SZ wrote and revised the manuscript. All authors gave feedback, read and approved the final manuscript.
The authors declare that they have no competing interests.
Consent for publication
Ethics approval and consent to participate
All animal work was approved by the University of Illinois at Urbana-Champaign Institutional Animal Care and Use Committee (IACUC), protocol numbers 15202 and 14209. Permits for bat collection and export were issued by the Republic of Trinidad and Tobago. No human subjects were used in this study.
Library quality control and note on bat aligment. (DOCX 29 kb)
Gene Expression Distributions Within & Between Species. Box and violin plots of FPKM (Fragments per kilobase per million reads) for each sample. 1e-3 FPKM was defined as the minimum detectable expression value. Violin plots (gray areas) of the density of reads have a bimodal distribution, one corresponding to genes with zero expression (below the cutoff) and the other to active genes. Boxplots indicate that all libraries have similar distribution of gene expression. The x-axis shows each individual sample used in the analysis. A: Forelimb and hindlimb of bats. B: Forelimb and hindlimb of mouse. C: Forelimb and hindlimb of pig. D: Forelimb and hindlimb of opossum. (TIF 7475 kb)
Analysis of consistency among replicates. Hierarchical clustering was used to determine similarity between replicates. A: All samples for bat. B: All samples for mouse. C: All samples for pig. D: All samples for opossum. (TIF 8736 kb)
Statistical significance of hierarchical clustering of all stages and limbs of each species (Fig. 2 in text). Red boxes are statistically significant clusters. A: Bat, B: Mouse, C: Opossum, D: Pig. The R package pvclust was used to determine statistical significance. FL = forelimb, HL = hindlimb. In A, St13, St14, and St15 correspond to the ridge, bud, and paddle stages respectively. In B, StW2, StW3-4, and StW6 correspond to the ridge, bud, and paddle stages. In C, St27 and St30, St28 and St31, and St29 and St32 correspond to the ridge, bud and paddle. For D, St20, St22, and St26 correspond to the ridge, bud and paddle stages. (TIF 1185 kb)
Statistical significance of hierarchical clustering of pairwise Spearmann coefficient values on fore (A-C) and hind (D-F) limbs (Fig. 3 in text). The pvclust R package was used to test for significance. No clusters were significant (p-value <=0.05). Abbreviations for stages are as in Supplementary Figure 3. (TIF 1270 kb)
Variation of transcriptome age index (VTAI) surrogate distribution for forelimbs. Bat, mouse, and pig forelimb trends are highly significant, while opossums are not. A: Bat B: Mouse C: Pig D: Opossum. (TIF 2530 kb)
Variation of transcriptome age index (VTAI) surrogate distribution for hindlimbs. As with the hindlimbs, bat, mouse, and pig forelimb trends are highly significant, while opossums are not. A: Bat B: Mouse C: Pig D: Opossum. (TIF 2474 kb)
Genes Divergently Expressed in Forelimbs vs. Hindlimbs. (DOCX 28 kb)
Genes with known roles in limb development (as identified by DAVID) that exhibit greatly divergent expression (75th percentile and above) among the limbs of all species. (DOCX 16 kb)
Hoxd12 WISH for opossum and bat forelimb and hindlimb. (TIF 23437 kb)
Additional replicates of Hoxd13 WISH for mouse, bat, and opossum forelimb and hindlimb. (TIF 28408 kb)
Additional replicates of Evx2 WISH for mouse and opossum forelimb and hindlimb. In mouse hindlimb ridges, ventral side is shown to highlight that purple staining representing Evx2 expression is not in the limb. (TIF 30642 kb)
Additional replicates of Hoxa13 WISH for mouse and opossum forelimb and hindlimb. (TIF 29860 kb)
Alignment of bat reads to the Myotis lucifugus genome using TOPHAT. FL and HL correspond to forelimb and hindlimb. (DOCX 29 kb)
Embryo samples used for whole-mount in situ hybridization. (DOCX 15 kb)
Accession Numbers & Primer Sequences For WISH. (DOCX 15 kb)
Genes with known roles in limb development that exhibit greatly divergent expression in the fore- and hind limbs of a single species and among the limbs of all species. (DOCX 20 kb)
Samples Used in Final Analysis. (DOCX 21 kb)
About this article
Cite this article
Maier, J.A., Rivas-Astroza, M., Deng, J. et al. Transcriptomic insights into the genetic basis of mammalian limb diversity. BMC Evol Biol 17, 86 (2017) doi:10.1186/s12862-017-0902-6
- Differential expression