这是indexloc提供的服务，不要输入任何密码
Skip to main content
BMC is moving to Springer Nature Link. Visit this journal in its new home.
Download PDF

Evolution of cetacean-specific conserved non-coding elements suggests their role in the limb changes during secondary aquatic adaptation

BMC Biology volume 23, Article number: 216 (2025) Cite this article

Abstract

Background

Limb morphology is particularly important for animals to inhabit different environments. Limb modifications (e.g., flipper-like forelimbs and hindlimb regression) are among the most critical secondary aquatic adaptation mechanisms enabling cetaceans to fully adapt to an aquatic environment. Exploring the molecular mechanisms underlying limb evolution in cetaceans has attracted considerable attention from evolutionary biologists.

Results

In the present study, conserved non-coding elements (CNEs) closely associated with limb development, which exhibited lineage-specific sequence divergence (nucleotide mutations and indels) in cetaceans, were identified using comparative genomics. These sequence divergences might have led to the loss of binding motifs for transcription factors involved in limb development and significant alterations in autoregulatory activity. A transgenic mouse was constructed to carry a cetacean-specific enhancer (i.e., hs1586), which exhibited a significant phenotypic difference in forelimb buds at embryonic day (E)10.5, supported by transcriptomic and epigenomic evidence. However, the phenotypic recovery after E11.5 suggested that enhancer redundancy in the mouse genome may have compensated for the effects caused by the incorporation of cetacean hs1586. This further suggests that the complex phenotypic changes of limbs in cetaceans are likely not driven by a single CNE but rather involve multiple CNEs and/or genes.

Conclusions

In summary, our study supports the functional role of CNE sequence divergence and the complex mechanisms underlying limb morphology changes in cetaceans.

Background

Under natural selection pressures, animals have evolved a wide variety of limb forms to adapt to different lifestyles (e.g., running, diving, digging, flying, climbing, swimming) [1, 2] and survival strategies (e.g., foraging, predation, predator avoidance, defense, reproduction) [3, 4]. Limb morphological changes in cetaceans (e.g., hindlimb reduction and flipper-like forelimbs) are critical evolutionary features linked to their secondary adaptation to an aquatic environment [5, 6]. Genetically, the signaling pathways (Wnt/β-catenin, Shh, and FGF) and gene regulatory mechanisms (e.g., Hox, Wnt, Fgf, and Tbx genes) that govern limb development in vertebrates exhibit remarkable consistency and high conservation [7]. Consequently, research on the morphological changes in cetacean limbs has long focused on functional genes [8]. It has been suggested that the loss of hindlimbs in the pantropical spotted dolphin (Stenella attenuata) is due to the disruption of Fgf8 expression and the absence of Hand2, yet the underlying causes of these abnormalities remain poorly understood [9]. Additionally, some genes involved in limb development, such as Hoxd11-Hoxd13 and LMBR1, have exhibited signs of accelerated evolution and sequence divergence in cetaceans [10,11,12], which were thought to potentially drive the limb morphological changes. However, due to the pleiotropy of gene functions [13,14,15], it remains uncertain whether the accelerated evolution of limb development-related genes in cetaceans has truly driven the changes pertaining to limb morphology.

Over the past decade, comparative genomics and bioinformatics studies have shown that non-coding regulatory elements (e.g., conserved non-coding elements, CNEs and enhancers) play a crucial role in limb development across diverse animal species [16,17,18,19]. During limb development, non-coding regulatory regions can bind to specific transcription factors (TFs) to regulate gene expression in a precise spatiotemporal manner [20, 21], thereby controlling cell fate and ensuring the accurate formation of limb morphology. Alterations in the sequences of non-coding regulatory regions could disrupt original spatiotemporal regulatory mechanisms [22, 23], providing new insights into the evolutionary development of limbs in different animals. In squamates (snakes and lizards), for instance, limb-related CNEs exhibit widespread sequence divergence [16, 18, 24]. Specifically, in extant snakes, a 17-bp specific deletion in the enhancer ZRS of the Shh has been identified as a key molecular mechanism underlying limb regression [25, 26]. However, the potential role of sequence divergence in non-coding regulatory elements in the morphological changes of cetacean limbs remains poorly understood.

In this study, a series of CNEs exhibiting specific sequence divergence (nucleotide mutations and indels) in cetaceans (hereinafter referred to as cetacean-specific CNEs) and their association with limb morphology changes were identified at the genome-wide scale. Studies have shown that CNEs are often candidate regions for regulatory elements, such as enhancers, promoters, and silencers [27, 28]. H3K27ac and H3K4me1 are histone modification marks associated with active enhancers, with H3K27ac specifically marking active regulatory elements to distinguish them from poised enhancers, while H3K4me1 is generally enriched in both active and poised enhancers [29]. We performed overlap analysis with public ChIP-seq data from the ENCODE database and identified cetacean-specific CNEs that may function as enhancers. Cellular functional experiments demonstrated that the enhancer activity of these cetacean-specific CNEs was significantly altered. Additionally, the transgenic mouse model carrying a cetacean-specific enhancer provided in vivo evidence supporting the impacts of sequence divergence in non-coding regulatory elements on limb development. This study could provide some new insights into the molecular mechanisms underlying cetacean limb changes in the secondary aquatic adaptation.

Results

Identification of CNEs associated with limb development in cetaceans through genome-wide screening

CNEs associated with limb development exhibit cetacean-specific sequence divergence

Thirty-eight representative mammals (26 marine and 12 terrestrial) were selected to screen CNEs (Fig. 1A and Additional file 1: Table S1). The high collinearity among the genomes of human (H. sapiens), mouse (M. musculus), and vaquita (P. sinus) ensured the quality of the identified CNEs (Additional file 2: Fig. S1a). A total of 333,341 CNEs (at least 50 bp in length) were screened, of which 3324 showed an accelerated evolutionary trend, and 2944 displayed specific indels in cetaceans (Additional file 2: Fig. S1b). The overlap analysis of the 6268 cetacean-specific sequence-divergent CNEs with the ChIP-seq data for H3K27ac and H3K4me1 modifications associated with limb development from the ENCODE database identified 745 and 1786 CNEs, respectively, which were most abundant during the limb bud initiation stage (E10.5) (Fig. 1B and Additional file 1: Tables S2, S3, and Additional file 2: Fig. S1c). Functional annotation of H3K27ac-modified CNEs (n = 745) using GREAT revealed that genes associated with these CNEs were significantly enriched in the Gene Ontology (GO) biological processes related to limb initiation and morphogenesis, including regulation of cartilage development (GO:0061035), embryonic limb morphogenesis (GO:0030326), limb development (GO:0060173), embryonic digit morphogenesis (GO:0042733), and anterior/posterior pattern specification (GO:0009952), among others. (Fig. 1C, and Additional file 1: Table S4, and Additional file 2: Fig. S2). GREAT functional annotation also revealed mammalian phenotype terms related to these CNEs from the Human and Mouse Gene Ontology Projects, with a focused on limb developmental abnormalities, including abnormality of finger (HP:0001167), abnormality of the ulna (HP:0002997), abnormal metacarpal bone morphology (MP:0003073), abnormal fibula morphology (MP:0002187), abnormal limb bud morphology (MP:0005650), and abnormal pelvic girdle bone morphology (MP:0004509) (Fig. 1D and Additional file 1: Table S4, and Additional file 2: Fig. S3). These results suggest that cetacean homologs of these CNEs may influence limb development. Additionally, GO and mammalian phenotype analyses indicated that these CNEs are associated with the development of the brain, heart, kidneys, and lungs. After filtering with public ChIP-seq data related to development of these organs from the ENCODE database (Additional file 1: Table S2), 163 CNEs with cetacean specific sequence divergence (hereinafter referred to as cetacean-specific CNEs) were eventually identified to be potentially related to limb changes of cetaceans (Additional file 1: Table S5 and Additional file 2: Fig. S4).

Fig. 1
figure 1

Identification of cetacean-specific CNEs closely associated with limb development. A A simplified phylogenetic tree of 38 representative species, with blue representing cetaceans (baleen whales and toothed whales). B The proportion of sequence-divergent CNEs in cetaceans associated with limb development, as reflected in H3K27ac-marked ChIP-seq data across different stages of mouse limb development. Blue represents indels, and orange represents accelerated evolution. C and D respectively display the GO biological processes and mammalian phenotype terms enriched in the CNEs with specific sequence divergence in cetaceans, as obtained through GREAT functional annotation. C Left: Bars depict Benjamini & Hochberg adjusted p-values derived by a one-sided Fisher’s exact test. Right: Many cetacean-diverged CNEs (dots) are far away from the transcription start site of genes in these sets. CNEs associated with the same gene have the same color. D Left: human phenotype terms. Right: mouse phenotype terms. Bars depict Benjamini & Hochberg adjusted p-values derived by a one-sided Fisher’s exact test

Full size image

Loss of motifs binding with transcription factors involved in limb development in cetacean-specific CNEs

To determine whether cetacean-specific CNEs alter the TFs they bind to, CNEs closely associated with limb development were selected as representatives. Firstly, TFs prediction was performed on CNE90, CNE227, CNE187, CNE239, CNE320, CNE315, CNE218, and CNE219, which are associated with the genes of paired-like homeodomain transcription factor 1 (Pitx1), bone morphogenetic protein 2 (Bmp2), transcription factor 4 (Tcf4), T-box transcription factor 15 (Tbx15), T-box transcription factor 3 (Tbx3), and Alx homeobox 4 (Alx4), respectively. These CNEs exhibited signals of accelerated evolution in cetaceans (Additional file 2: Figs. S5, S6). The potential TFs that bind to these CNEs were predicted using two databases, CIS-BP and AnimalTFDB 4.0 (Additional file 1: Tables S6, S7). Subsequently, based on the Spearman correlation analysis, TFs significantly correlated with the expression of potential target genes were identified (Fig. 2A and Additional file 1: Table S8). GO enrichment analysis of these TFs was performed using Metascape, which showed that biological processes significantly enriched in the TFs predicted by mouse CNEs (including forelimb morphogenesis, cartilage development, and regulation of osteoblast differentiation) were not found in the GO enrichment analysis of TFs predicted by cetacean CNEs (Fig. 2B and Additional file 1: Table S9). Secondly, CNE622 (hs1262), CNE682, CNE497, CNE531, and CNE629, which are associated with the genes of short stature homeobox 2 (Shox2), homeobox gene A13 (HoxA13), paired box gene 9 (Pax9), Wnt family member 5A (Wnt5A), and bone morphogenetic protein receptor 1B (Bmpr1B), respectively, have fragment deletions specific to cetaceans (Additional file 2: Fig. S7). Predictive analysis based on the JASPAR database reveals that key TFs (e.g., Pitx1, Twist2, Myod1, Sox10) regulating limb development are likely to bind at the deletion sites of those cetacean homologous CNEs (Fig. 2C and Additional file 1: Table S10). Additionally, the enhancer hs1586 located in the intron of the Gli3, which exhibited a 17-bp deletion and an accelerated evolutionary signal in cetaceans (Additional file 2: Fig S8), was predicted to be bound by key TFs regulating limb development, such as Shox2, Hoxd9, Ets1, and Pitx1, at its deletion site (Fig. 2D and Additional file 1: Table S10). Taken together, the widespread sequence divergences in cetacean CNEs may lead to the loss of binding sites for TFs regulating limb development.

Fig. 2
figure 2

Prediction of TFs binding for cetacean-specific CNEs and mouse CNEs. The CNEs in cetaceans have undergone sequence divergence, which may lead to changes in the TFs that bind to them. To validate this hypothesis, we used mouse CNEs as controls and randomly selected cetacean CNEs to perform TFs predictions for both cetacean and homologous mouse CNEs. These CNEs are associated with key genes involved in limb development. A CNEs associated with Bmp2, Pitx1, Tbx15, Tbx3, and Tcf4 were used to predict potential TFs. Spearman correlation analysis was then applied to identify TFs that were significantly co-expressed with their target genes. Significantly negatively correlated: − 1 < Spearman’s rank correlation coefficient <  − 0.5. Significantly positively correlated: 0.5 < Spearman’s rank correlation coefficient < 1. P-values were adjusted using the Holm-Bonferroni method. B GO functional enrichment analysis using Metascape on TFs that are significantly co-expressed with their target genes. (top) TFs predicted for mouse CNEs, (bottom) TFs predicted for dolphin CNEs. Bars depict Benjamini & Hochberg adjusted p-values derived by a one-sided Fisher’s exact test. C, D CNE629, CNE682, CNE622, CNE531, and hs1586 are associated with Bmpr1b, Hoxa13, Shox2, Wnt5a, Pax9, and Gli3, respectively, and exhibited cetacean-specific indels. Based on the JASPAR database, potential TFs that may bind to the indel sites of these CNEs in cetaceans were predicted. C Top: Bmpr1b-CNE629; middle: Hoxa13-CNE682, Shox2-CNE622; bottom: Wnt5a-CNE531, Pax9-CNE497. D Gli3-hs1586

Full size image

Cetacean-specific sequence divergence affects the regulatory activity of CNEs

Functional annotation of cetacean-specific CNEs using GREAT revealed associated genes (e.g., Gli3, Pax9, HoxA13, Alx4, Tcf4, Bmp2, Pitx1, Shox2), which play important roles in the limb bud initiation and limb development (Fig. 3A and Additional file 1: Table S5). Dual luciferase assays for target CNEs and enhancers showed that the regulatory activities of cetacean homologous sequences of hs1586, hs1262, CNE682, CNE497, CNE218, CNE219, CNE187, and CNE227 were significantly reduced compared to those of mouse (Fig. 3B–D and Additional file 2: Figs. S9a–S9e). However, compared to the mouse, CNE90 exhibited much higher regulatory activity in the minke whale and much lower regulatory activity in the bottlenose dolphin (Additional file 2: Fig. S9f), which is consistent with the difference in nucleotide substitutions of CNE90 between the two species. In addition, compared to the mouse, other limb development-related CNEs differentiated in cetaceans (e.g., CNE320, CNE315, CNE239, CNE531, CNE629) exhibited significant differences in regulatory activities (Additional file 2: Figs. S9g–S9k).

Fig. 3
figure 3

The regulatory activity of cetacean-specific CNEs was evaluated using dual-luciferase assays. A Simplified gene regulatory signaling pathways for early limb development in vertebrates, with genes associated with cetacean-specific CNEs playing important roles in the initiation and growth of limb development. B–D Compared with the positive control mouse, the dual-luciferase assay sequentially validated the regulatory activity changes of hs1586, CNE622, and CNE682 in cetaceans. From left to right, the groups are Control, M. musculus, B. acutorostrata, and T. truncatus. B hs1586, C CNE622, D CNE682. E The dual-luciferase assay validated changes in the regulatory activity of mouse hs1586 carrying cetacean mutations. From left to right, the groups are Control, M. musculus, T. truncatus, and the mutant group. F–H The dual-luciferase assay validated changes in the regulatory activity of mouse hs1586 and mouse hs1586 carrying cetacean mutations after transcription factor binding. From left to right, the groups are: M. musculus without transcription factor, M. musculus with transcription factor, and the mutant with transcription factor group. F Ets1, G Hoxd9, H NR4A2. Data are meaniption factor, and the mutant with transcription P-value was calculated using Studentnd t-test (*, P < 0.05; **, P < 0.01; ***, P < 0.001)

Full size image

To further test whether sequence divergence in cetacean CNEs leads to the loss of TF binding sites and thus affects autoregulatory activity, hs1586 was used as an example for experimental validation. Initially, a dual-luciferase assay was performed using the plasmid carrying mouse hs1586 and the cetacean mutant hs1586. The results showed that the regulatory activity of hs1586 was significantly reduced after mutation (Fig. 3E), suggesting that the deletion of a fragment in the cetacean hs1586 significantly affected its regulatory activity. Subsequently, the dual-luciferase assays were performed using plasmid containing TFs (NR4A2, Ets1, Hoxd9) with the plasmid carrying the mouse hs1586 and the cetacean mutant hs1586, respectively. The results showed that the regulatory activities of both mouse and cetacean mutant hs1586 were significantly upregulated after adding TFs. However, the regulatory activity of the mutant hs1586 remained significantly lower than that of the mouse hs1586 (Fig. 3F–H). These results suggested that the specific fragment deletion of cetacean hs1586 does lead to reduced regulatory activity due to its inability to bind with corresponding TFs.

In vivo evidence supporting the role of cetacean-specific hs1586 in forelimb bud morphology

Based on the cetacean-specific sequence divergence of hs1586 (including an accelerated evolutionary signal and a 17-bp deletion) and its significantly reduced enhancer activity, we knocked the homologous cetacean hs1586 sequence into the mouse (C57BL/6 J) genome to generate a transgenic mouse model (Fig. 4A and Additional file 2: Fig. S10). Phenotypic analysis revealed that at E10.5 only, the forelimb buds of transgenic homozygous mice (THM) (right: mean = 2.14882 mm2, left: mean = 2.19051 mm2) were significantly larger than those of wild-type mice (WT) (right: mean = 1.99005 mm2, left: mean = 2.00774 mm2) (right p-value = 0.0478, left p-value = 0.0270) (Fig. 4B and Additional file 1: Tables S11–S12 and Additional file 2: Fig. S11). However, no significant phenotypic differences were observed at other developmental stages (i.e., E11.5, E18.5, and postnatal day 1 (P1)) (Additional file 1: Tables S11–S16 and Additional file 2: Figs. S11, S12, S13). RNA-seq analysis showed that, compared to WT, 674 and 159 significantly differentially expressed genes (DEGs) were detected in forelimb and hindlimb buds at E10.5, respectively, which was much more than those at E11.5 and E12.5 (Fig. 4C and Additional file 1: Tables S17–S19 and Additional file 2: Fig. S14). At E10.5, the DEGs in forelimb buds included many key genes that regulate limb development, such as Shh, Fgf4, Fgf8, HoxA13, and HoxD13 (Fig. 4D, and Additional file 1: Table S17), with only a few DEGs involved in limb development in hindlimb buds (Additional file 1: Table S17 and Additional file 2: Fig. S15a). Additionally, GO and KEGG enrichment analyses revealed that DEGs in forelimb buds at E10.5 were significantly enriched in biological processes and signaling pathways closely associated with limb development (Figs. 4E, F). In contrast, among the signaling pathways enriched in the DEGs of the hindlimb buds, only the TGF-β, MAPK, and Wnt signaling pathways were significantly associated with limb development (Additional file 2: Figs. S15b, S15c). RT-qPCR further confirmed upregulation of Sox6, Tbx1, Dlx5, HoxD13, Runx3, Runx1, and Pax9, and downregulation of Pax3, Fgf8, and Fgf4 in the forelimb buds at E10.5 compared with WT (Additional file 2: Fig. S15d).

Fig. 4
figure 4

In vivo molecular evidence supports the morphological changes in the forelimb buds of THM at E10.5. A A transgenic mouse model was generated by knocking in the homologous cetacean hs1586. B The comparison of limb bud size between WT and THM. Data are mean standard error of three independent measurements. P-value was calculated using Student oft-test (*, P < 0.05; **, P < 0.01). top: forelimb buds, bottom: hindlimb buds. Scale bar, 1:2500 µm. C The number of DEGs between WT and THM during E10.5–E12.5. D A volcano plot illustrating differentially regulated gene expression from RNA-seq analysis in forelimb buds between WT and THM at E10.5. Genes that are upregulated and downregulated are shown in red and blue, respectively. Values are presented as the log2 of tag counts. E GO functional clustering of genes that were significantly altered in forelimb buds of THM for biological processes. F KEGG pathway analysis of DEGs in forelimb buds of THM transcriptome

Full size image

Epigenetic analysis revealed that within 250 kb upstream and downstream of hs1586 (chr13:15,548,974–15,550,831 mm10) at E10.5, WT had 31 activated enhancers, while THM had only 7 (Fig. 5A and Additional file 1: Table S20). The scanned activated enhancers were mapped to their respective genomes, and it was found that 17 activated enhancers were located in the intronic region of the Gli3 in WT, consistent with previous results (Additional file 1: Table S21). However, only 4 activated enhancers were localized to the Gli3 intronic regions in THM (Fig. 5B). Subsequently, at E11.5, epigenetic analysis revealed that WT had 31 activated enhancers, 21 of which were located in the Gli3 intron region. For THM, there were 23 activated enhancers, and 12 of which were located in the Gli3 intronic region (Fig. 5C, D, and Additional file 1: Table S20). Interestingly, there was no significant difference in the transcription level of the predicted target gene Gli3 (log2FC =  − 0.5168, FDR = 0.1277) compared with the WT at E10.5, while the transcription level of Inhibin beta A (Inhba) (log2FC = 4.2430, FDR = 8.77E-11), a member of the transforming growth factor-β family that was located 461 kb downstream of the target region (chr13:15,548,974–15,550,831 mm10), was significantly upregulated (Additional file 2: Fig. S16a). Notably, Inhba expression was not detected during hindlimb bud development at E10.5. Moreover, at E11.5 and E12.5, Inhba expression in both forelimb and hindlimb buds showed no significant difference compared with WT (Additional file 2: Figs. S16a, S16b), suggesting that the absence of Inhba expression in the hindlimb bud at E10.5 may be related to normal hindlimb bud development.

Fig. 5
figure 5

Epigenetic analysis of forelimb buds at E10.5 and E11.5 in THM and WT. A Normalized H3K27ac (red) and H3K4me1 (blue) epigenetic signals in the region spanning 250 kb upstream and downstream of the hs1586 locus for mixed samples from each group in the forelimb buds of THM and WT at E10.5. Top: WT, bottom: THM. B At E10.5, activated enhancers within 250 kb upstream and downstream of the hs1586 in WT and THM were mapped to the genome, respectively. Top: WT, bottom: THM. C Normalized H3K27ac (red) and H3K4me1 (blue) epigenetic signals in the region spanning 250 kb upstream and downstream of the hs1586 locus for mixed samples from each group in the forelimb buds of THM and WT at E11.5. Top: WT, bottom: THM. D At E11.5, activated enhancers within 250 kb upstream and downstream of the hs1586 in WT and THM were mapped to the genome, respectively. Top: WT, bottom: THM

Full size image

Discussion

Changes in regulatory activity of cetacean-specific CNEs associated with limb development are the mechanisms underlying limb changes of cetaceans

Sequence divergence of non-coding regulatory elements has been regarded as playing an important role in the evolution of vertebrate limb morphology [16, 18, 30]. However, it remains unclear whether this holds true in cetaceans and how the evolution of CNEs could contribute to the limb morphological changes in cetaceans. In this study, we found that CNEs associated with cetacean limb development have undergone species-specific sequence divergence, as revealed by bioinformatics analyses. This pattern is similar to that observed in other species, such as snakes and certain lizards, with drastic limb morphological changes [16, 18, 24, 31]. Our analysis showed that cetacean-specific CNEs associated with limb development account for the highest proportion during the early limb bud stage (E10.5), suggesting their potential role in early limb development. A previous study has shown that the regression of the hindlimbs in bottlenose dolphins is associated with the disruption of Fgf8 expression, leading to the failure of the apical ectodermal ridge (AER) maintenance and the inability to establish the zone of polarizing activity (ZPA) [9]. Similarly, in snakes and lizards, limb reduction is often accompanied by early necrosis of the AER [26, 32, 33]. We speculate that cetacean limb morphological changes may likewise occur during early developmental stages. TFs prediction analysis determined that these cetacean-specific CNEs have lost the binding motifs of key TFs (e.g., Pitx1, Shox2, Twist2, Ets1, and Sox10) involved in limb development [34,35,36]. This may disrupt the original TFs binding and alter the precise spatiotemporal regulation of relevant genes, as demonstrated in previous studies [37]. However, compared to snakes and lizards, cetacean-specific CNEs exhibit different patterns of sequence divergence, which may lead to specific reshaping of their regulatory networks during evolution. Whether this change specifically drives the evolution of cetacean limb morphology remains to be further investigated.

The influence of CNEs on the morphological development of tissues and organs stems from their precise spatial and temporal regulation of genes [38,39,40,41]. Point mutations and indels within CNEs can alter their cis-regulatory activity and then lead to the change of their regulatory function [24, 42, 43], thereby affecting the development of tissues and organs. Our analysis revealed that genes associated with cetacean-specific CNEs play important roles in limb initiation and limb bud formation (Fig. 3A). For example, abnormal expression of Gli3 can lead to several human genetic disorders involving polydactyly and syndactyly, such as Greig cephalopolysyndactyly syndrome (GCPS), Pallister-Killian syndrome (PHS), and postaxial polydactyly [44]. Additionally, in a mouse model, knockout of Alx4 leads to polydactyly [45], while knockout of Shox2 results in humeral shortening [46]. Therefore, significant changes in the regulatory activity of cetacean-specific CNEs may alter the expression patterns or expression levels of their potential target genes, thereby affecting limb morphology during development. The functional significance of cetacean-specific CNEs was further supported by the transgenic mouse model experiments. Compared with WT, THM that carried the cetacean-specific hs1586 exhibited significant phenotypic differences in the forelimb buds at E10.5, associated with significant in vivo transcriptomic differences (including DEGs at the same developmental stage and their corresponding GO/KEGG enrichment analyses). However, unexpectedly, the transcriptional expression level of Gli3 did not show a significant change in THM, which may be due to the presence of multiple regulatory CNEs and enhancers within its introns [47, 48]. A similar situation is observed for the Hox and Fgf8 genes, which are finely regulated by multiple cis-regulatory elements during limb development [49, 50]. This regulatory redundancy may provide a buffering effect to some extent. Future studies combining Hi-C and ChIP-seq could help elucidate the impact of active CNEs and enhancers within the Gli3 introns on its expression.

Thus, the loss of potential TFs binding and the significant changes in regulatory activity caused by the cetacean-specific CNEs may be the mechanism for the limb morphology evolution in cetaceans.

The phenotypic recovery in the transgenic mice suggested a complex mechanism for the limb changes of cetaceans

Although the significant phenotypic differences in the forelimb buds of THM at E10.5 suggest that cetacean-specific mutated CNEs may be the molecular mechanisms underlying cetacean limb morphological changes, this phenotypic difference was recovered after E11.5. Epigenetic analysis of the forelimb buds showed that, compared to E10.5, the number of potential activated enhancers near hs1586 in THM increased at E11.5, suggesting the presence of redundant CNEs in mouse genome that are possibly absent in cetacean genome. Like the compensatory mechanism (i.e., functional redundancy among enhancers) revealed in previous studies [46, 49,50,51], the present finding of redundant CNEs in the mouse genome could, to some degree, provided explanation for the phenotypic recovery in the later limb development in the THM. In other words, the morphological changes in cetacean limbs, as a complex trait, might involve multiple genes or regulatory elements and therefore cannot be revealed through evolutionary analyses and functional experiments that focused solely on a single or very few CNEs or enhancers. Vertebrate limb development is governed by a tightly regulated and complex gene regulatory network composed of multiple signaling pathways, such as Shh, Fgf, Bmp, and Wnt [52,53,54]. Our analysis revealed that genes associated with cetacean-specific CNEs play roles in several key pathways. Gli3 and Alx4 regulate limb bud initiation and patterning by participating in the Shh signaling pathway [45, 55]; BMPs and the Grem1 signaling pathway control the proliferation and apoptosis of interdigital cells, which is closely related to the formation of webbed or fin-like limbs [8]; and Wnt5a, through the Wnt signaling pathway, interacts with the Fgf signaling pathway to coordinate mesenchymal cell arrangement and limb bud patterning during early development [56]. These findings suggest that cetacean-specific CNEs may influence limb developmental patterns and morphology in cetaceans by regulating the expression of genes in key signaling pathways. When a gene is located in a densely regulated region, its expression is typically controlled by multiple enhancers acting in concert, further highlighting the complexity of limb developmental regulatory networks [48,49,50]. Therefore, future studies could adopt a stepwise approach to edit multiple cetacean-specific enhancers or CNEs associated with the same gene, or to progressively target cetacean-specific enhancers or CNEs across multiple genes, in order to construct more complex transgenic mouse models. This would enable a deeper investigation into how sequence divergence in non-coding regulatory elements has contributed to the evolution of cetacean limb morphology.

Conclusions

This study demonstrated that non-coding regulatory elements in cetaceans with specific sequence divergence may play a role in influencing the limb morphology during development, which could provide some novel insights regarding the molecular mechanisms underlying the evolution of cetacean limb changes. However, the limited availability of high-quality cetacean genome and transcriptome data remains a major limitation of this study. In the future, global collaborations, increased sequencing samples, and the integration of systematic bioinformatics, comparative genomics, and Evo-Devo evidence are essential for further in-depth understanding of the molecular regulatory mechanisms driving cetacean limb evolution.

Methods

Identification of CNEs

Firstly, prior to the multi-genome alignment of 38 representative species, a simple collinearity analysis of the human (H. sapiens), mouse (M. musculus), and vaquita (P. sinus) genomes was performed using TBtools software [57]. This analysis assessed the collinearity between the cetacean and human or mouse genomes, ensuring the high quality of the screened CNEs. Subsequently, based on the Latin names of 38 representative species, a phylogenetic tree was acquired from the TimeTree website (http://www.timetree.org/) [58]. Whole-genome multiple alignments for the 38 species were obtained from the LASTZ v1.04.03 (parameters: H = 2000, K = 2400, L = 3000, Y = 9400) pairwise alignments [59], with the human genome (GRCh38/hg38) as reference. The UCSC tools were applied for “chaining” and “netting” to ensure that all the alignments were non-overlapping. To ensure accurate sequence alignment, we referred to the parameters described in the genomic alignment analysis by Hecker and Hiller [60]. Then, Multiz (version 11.2) was used to build a multiple genome alignment using human as the reference [61]. PhyloFit of the PHAST v1.5 software package was used to estimate a neutral model using fourfold degenerate (4d) sites that were extracted from the whole-genome alignments [62]. PhastCons (parameters: target-coverage = 0.3, expected-length = 45, rho = 0.3) was then performed to identify the evolutionarily conserved elements [62]. Raw conserved elements were merged when their distance was less than 10 bp. Finally, CNEs of at least 50 bp were retained for subsequent analysis.

Identification of cetacean-specific CNEs

To identify CNEs that exhibit accelerated evolution in cetacean-specific branches, we detected shifts in the DNA substitution rate of genomic regions using PhyloAcc software and identified cetacean-specific accelerated genomic elements from a set of CNEs [63]. We compared the marginal likelihood under three models: a null model (M0), where all branches are neutral or conserved; an accelerated model (M1), where branches leading to cetaceans are accelerated; and a full model (M2) with no constraints on latent states. We defined a conserved state as having a substitution rate lower than the background branch lengths (r1 < 1), and an accelerated state as having a substitution rate higher than the conserved state (r2 > r1). Two Bayes factors (BF) were used: BF1 (comparing M1 with M0) and BF2 (comparing M1 with M2) to identify elements that are accelerated exclusively in cetacean lineages. BF1 captures elements accelerated in cetacean groups regardless of the rest of the clade’s evolutionary pattern, while BF2 is crucial for excluding elements that show acceleration in other species. Cetacean-accelerated elements were defined as those with BF1 ≥ 5 and BF2 ≥ 2. Additionally, we used in-house Perl scripts to identify the CNEs with species-specific indels in cetaceans, which however are highly conserved within other mammals. This is based on previous studies demonstrating that mutations (including substitutions and Indels) within CNEs may be associated with the phenotypic diversity of species, such as the loss of limbs in snakes [25].

Function annotation of CNEs

Temporal ChIP-seq data for H3K27ac and H3K4me1 modifications across different stages of mouse limb development were obtained from the ENCODE database. Using the “Bedtools” software (version 2.31.1), we performed an overlap analysis between the coordinates of mouse orthologous sequences of cetacean-specific CNEs and the publicly available ChIP-seq data, ultimately identifying cetacean-specific CNEs closely associated with limb development.

The genomic coordinates of CNEs were extracted using the human genome (GRCh38/hg38) as the reference, and the “Basal plus extension” option (up to 1,000 kb) in GREAT (version 4.0.4 http://great.stanford.edu) was applied to associate CNEs with nearby genes [64]. Then, we compared genes associated with cetacean- or control-accelerated CNEs to genes near all CNEs (background) and searched for any functional enrichment in GO biological processes and mammalian phenotypes from the Human Gene Ontology Project. To ensure the accuracy of functional enrichment analysis, GREAT utilizes two methods for multiple hypothesis testing correction: Bonferroni correction and False Discovery Rate (FDR) correction. Further details can be found on the GREAT website. We only retained annotation terms that contained more than five genes in total, including at least two genes associated with accelerated CNEs and at least 1.5-fold enrichment of tested CNEs over all CNEs.

TFs prediction and correlation analysis

Using AnimalTFDB-4.0 [65] and CIS-BP (threshold ≥ 8) [66], potential TFs that may bind to selected mouse CNEs and their orthologous cetacean-specific CNEs were predicted, with mouse as the reference species. To ensure prediction accuracy, only TFs that were consistently identified by both databases were retained. Subsequently, we collected transcriptomic data of mouse limbs at different developmental stages from the ENCODE database and analyzed the co-expression patterns between the predicted TFs and target genes using Spearman’s correlation [67], in order to determine whether these TFs are significantly correlated with the potential target genes regulated by CNEs. Additionally, we conducted functional enrichment analysis of the significantly co-expressed TFs using Metascape. For CNEs that are highly conserved in other mammals but exhibit cetacean-specific indels, we performed transcription factor predictions at their indels sites, utilizing the JASPAR public database [68].

Dual-luciferase assays

CNEs fragments were chemically synthesized and cloned into the pGl3.0-sv40-luciferase vector. HEK293T cells were cultured in a Dulbecco’s Modified Eagle Medium (DMEM) containing 10% fetal bovine serum at 37°(DMEM) c2, and 100% humidity. Cells were seeded into 24-well plates. CNE constructs were transfected 24 h after seeding, using Polyethylenimine Linear (PEI) MW25000(YEASEN, #40815ES08) as the transfection reagent according to the manufacturer’s instructions. Two microliters of transfection reagent and 500 ng of the constructed pGl3.0-sv40-luciferase were mixed in 50 µL of basic medium. After sitting for 5 min at room temperature (RT), the above mixture was co-transfected with 50 ng of phRL-sv40-Rluc and then co-transfected into HEK293T cells after an incubation period of 20 min at RT. Cellular lysates were collected 48 h after transfection using the Dual Luciferase Reporter Assay Kit (Vazyme, #DL101). Luciferase activity was measured using the Dual-Luciferase Reporter Assay System (Promega) on a multi-mode microplate reader. The assay for each CNE was repeated three times. Luciferase activity was calculated as the ratio of firefly luciferase activity to Renilla luciferase activity in each well. Significance was assessed by Student’s t-test.

Generation and identification of the transgenic mouse model

We knocked the cetacean-specific hs1586 into mouse via CRISPR/Cas9 system. Firstly, two sgRNAs (gRNA1: 5′ CCAATATGCTTTATCTATGC 3′, PAM: AGG, gRNA2: 5′ AGTCAGTGAAGCACTGTAAT 3′, PAM: CGG, chemical synthesis) targeting sequences near the insertion site were designed and transcribed in vitro, and the donor vector with the inserted fragment was designed and constructed in vitro. Secondly, Cas9 mRNA (GemPharmatech Co., Ltd), sgRNA, and donor vector were co-injected into zygotes. Thereafter, the zygotes were transferred into the oviduct of pseudopregnant ICR females at 0.5 days post-coitum. All the offsprings of ICR females (F0 mice) were identified by PCR and crossing positive F0 mice with WT to build up heterozygous mice. Finally, specific primers were designed based on the cetacean hs1586 and mouse hs1586 sequences to genotype the offspring of heterozygous knock-in mice and identify homozygous knock-in mice. Cetacean hs1586 primers: Forward 5′GCGCTTACTTTTGAGAAAAGATCGAGTCCAGA3′, Reverse 5′GCACAGATAGTCACACCACTGCACATAGTTGA3′. Mice hs1586 primers: Forward 5′ATGAGCACTTCTTTCTATTGGCTGATGGTAGTGGA3′, Reverse 5′GTTTAGTAAGACATCACGTCATCCACTCAACTGCA3′.

All mice in this study were raised in the SPF animal room of the Animal Resources Center of Nanjing Normal University (NJNU-ARC), and the management of experimental animals was conducted under the guidance of the Animal Use and Ethics Committee of NJNU (ethics approval number: IACUC-20200501). Mice were maintained under a standard 12-h light/dark cycle and monitoring according to NJNU-ARC policies and procedures. During sample collection, we strictly followed the guidelines issued by the International Council for Laboratory Animal Science (ICLAS) for the euthanasia of mice. The mice were placed in a sealed chamber with a gradual introduction of carbon dioxide (CO2) to induce euthanasia. The absence of respiration and other vital signs, such as chest movement, was observed to confirm death, after which CO2 flow was maintained for an additional 2 min to ensure complete euthanasia.

Phenotypic analysis

For embryos at stages E10.5, E11.5, and E18.5, we strictly followed the mating time of the female and male mice for sampling. Specifically, the female and male mice were housed together at 5 PM on the first day, and vaginal plugs were checked before 9 AM the next day, which was designated as E0.5. Embryos were then collected at 6 PM on days 10, 11, and 18 of pregnancy, and these embryos were considered to be at developmental stages E10.5, E11.5, and E18.5, respectively. The day of mouse birth is designated as postnatal day 0 (P0), and the following day is referred to as postnatal day 1 (P1). All sampling was conducted according to the aforementioned time points in this study. Limb buds at E10.5–E11.5 were separated along the body axis under Leica Microsystems (M165FC) to remove as much excess muscle tissue as possible. The isolated limb buds were photographed at the same magnification (Leica Microsystems Ltd M165FC 1 ×). To characterize the limb buds for different samples, the area of limb buds was measured using ImageJ software (http://rsbweb.nih.gov/ij/) at the same scale bar (1:2500 µm) [69]. For embryos at E18.5 and postnatal day one, alcian blue/alizarin red (Sangon Biotech (Shanghai) Co., Ltd. A600298, A506786) was used to stain the cartilage and osteogenesis of the limbs [70]. Again, pictures were taken at the same magnification (Leica Microsystems Ltd M165FC 0.73 ×), and the same scale bar was used to measure osteogenesis length. All the above measurements were performed in nine independent replicates. To facilitate the t-test (the significance is indicated as *p < 0.05), each set of measurements was tested for normal distribution, and the results are presented in the supplementary material (Additional file 1: Tables S11–S16).

RNA sequencing analysis

E10.5–E12.5 embryos from WT and THM were collected from different pregnant mice and divided into three independent replicates. In each replicate, the same number of limb buds was pooled from littermate embryos. All procedures related to RNA extraction, library construction, sequencing, and transcriptome data analysis (including quality control, read alignment, differential expression analysis, and functional annotation) were performed by Beijing Biomarker Technologies Co., Ltd (http://www.biomarker.com.cn). Briefly, total RNA of each sample was extracted according to the instruction manual of the TRIzol Reagent (Life technologies, California, USA). RNA integrity and concentration were checked using an Agilent 2100 Bioanalyzer (Agilent Technologies, Inc., Santa Clara, CA, USA). The mRNA was isolated by NEBNext Poly (A) mRNA Magnetic Isolation Module (NEB, E7490). The cDNA library was constructed following the manufacturer’s instructions of NEBNext Ultra RNA Library Prep Kit for Illumina (NEB, E7530) and NEBNext Multiplex Oligos for Illumina (NEB, E7500). The constructed cDNA libraries were sequenced on a flow cell using an Illumina HiSeq™ sequencing platform. The raw sequences were processed using Fastp software (parameters: -i L006_R1_001.fq.gz -o L006_good_1.fq -I L006_R2_001.fq.gz -O L006_good_2.fq -Q -y -g -Y 10 –adapter_fasta adapter.fa -l 100 -b 150 -B 150 -j L006.json -h L006.html -w 6) to remove low-quality reads, including those containing only adapter sequences, those with unknown nucleotides (> 5%), or those with a Q20 score < 20% (indicating a sequencing error rate > 1%). The remaining high-quality reads were retained as clean reads. Subsequently, the clean reads were mapped to the mouse genome (mm10, GRCm38) with high precision using Hisat2 (version 2.0.4, parameters: –dta -p 6 –max-intronlen 5,000,000) to obtain their genomic alignment information. Gene expression levels were quantified using StringTie (version 2.2.1, parameters: –merge -F 0.1 -T 0.1) and normalized using FPKM (fragments per kilobase of exon per million mapped reads) as a metric for assessing transcript or gene expression levels [71]. DESeq [72] and q-value were employed and used to evaluate differential gene expression between WT and transgenic mice. After that, gene abundance differences between those samples were calculated based on the ratio of the FPKM values. The false discovery rate (FDR) control method was used to identify the threshold of the p-value in multiple tests in order to compute the significance of the differences. Here, only genes with an absolute value of log2 ratio ≥ 2 and FDR significance score < 0.05 were used for subsequent analysis. Genes were compared against various protein databases by BLASTX, including the National Center for Biotechnology Information (NCBI) non-redundant protein (Nr) database, Swiss-Prot database with a cut-off e-value of 1E-5. Furthermore, genes were searched against the NCBI non-redundant nucleotide sequence (Nt) database using BLASTn (parameters: -nohead -evalue 1E-5 -tophit 1 -m 0 -topmatch 1) with a cut-off e-value of 1E-5. Genes were retrieved based on the best BLAST hit (highest score) along with their protein functional annotation. To annotate the gene with gene ontology (GO) terms, the Nr BLAST results were imported into the Blast2GO program [73]. This analysis mapped all of the annotated genes to GO terms in the database and counted the number of genes associated with each term. A Perl script was then used to plot GO functional classification for the unigenes with a GO term hit to view the distribution of gene functions. The obtained annotation was enriched and refined using TopGo (R package). The gene sequences were also aligned to the Clusters of Orthologous Group (COG) database to predict and classify functions [74]. KEGG pathways were assigned to the assembled sequences using the R package clusterProfiler (version 4.4.4) with the following parameters: minGSSize = 1, maxGSSize = 10,000, pAdjustMethod = “fdr”.

Quantitative reverse transcription PCR (RT-qPCR)

E10.5 embryos from WT and THM were collected from different pregnant mice. Forelimb buds were isolated along the spine and used for RT-qPCR. Each sample group consisted of an equal number of forelimb buds pooled from littermate embryos. The experiment was performed with three independent biological replicates. Total RNA was extracted according to the instruction manual of the RNA isolater Total RNA Extraction Reagent (Vazyme R401-01) and reverse-transcribed into cDNA using HiScript II Q RT SuperMix for qPCR (+ gDNA wiper) (Vazyme R223-01) according to the manufacturer’s instructions. Expression of mRNAs was determined using ChamQ SYBR qPCR Master Mix (2 × without ROX) (Vazyme Q321-02). The relative expression level of the target was calculated using the comparative Ct method. GAPDH was used as an internal control to normalize sample differences. Significance was assessed by t-test (*p < 0.05, **p < 0.01). The sequences of the primers used for RT-qPCR analysis are presented in Additional file 1: Table S22.

Epigenetic analysis (CUT&Tag)

E10.5–E11.5 embryos from WT and THM were collected from different pregnant mice. For each sample used for sequencing, an equal number of forelimb buds from different pregnant mice offspring was pooled. The genotype of all embryos was recharacterized prior to the experiments. H3K27ac and H3K4me1 [29, 75], which are closely related to enhancer activity, were selected to perform CUT&Tag. The overlapping regions of H3K27ac and H3K4me1 signal peaks were used to define activated enhancer.

CUT&Tag assays were performed by Wuhan IGENEBOOK Biotechnology Co., Ltd (http://www.igenebook.com) [76]. Briefly, 50,000 cells were harvested and incubated with concanavalin A coated-magnetic beads for 15 min at RT. Then, bead-bound cells were resuspended and incubated with the appropriate primary antibody (H3K4me1, CST5326, H3K27ac, ab4729) and IGG (2729 s, CST) overnight at 4 ℃. Afterwards, Goat anti-rabbit secondary antibody (ab206, Abcam) was added to this mixture for 30 min at RT. The magnet stand was used to remove unbound antibodies. pA-Tn5 adapter complex was added at this step at RT for 1 h, and then the unbound pA-Tn5 protein was removed. Next, cells were resuspended in tagmentation buffer at 37 ℃ for 1 h. After that, immunoprecipitated DNA was saved. To amplify libraries, 21 µL DNA was mixed with 2 µL of a universal i5 and a uniquely barcoded i7 primer, using a different barcode for each sample. A volume of 25 µL NEB Next HiFi 2 × PCR Master mix was added and mixed. The sample was placed in a Thermocycler with a heated lid using the following cycling conditions: 72 °C for 5 min (gap filling); 98 °C for 30 s; 14 cycles of 98 °C for 10 s and 63 °C for 30 s; final extension at 72 °C for 1 min and hold at 8 °C. Library was sequenced on Illumina NovaSeq 6000 with PE 150 method. Trimmomatic (version 0.36) was used to filter out low-quality reads [77]. WT clean reads were mapped to the mm10 genome by BWA (version 0.7.15), while THM clean reads were mapped to the mm10 genome-1 which had undergone sequence substitution in the region of chr13:15,548,974–15,550,831 [78]. Samtools (version 1.3.1) was used to remove potential PCR duplicates [79]. MACS2 software (version 2.1.1.20160309) [80] was used to call peaks by default parameters (bandwidth, 300 bp; model fold, 5, 50; p-value, 0.00001). Based on the analysis principle of Deseq [81], the enrichment ratio of H3K27ac and H3K4me1 was calculated in the genomes of WT and THM.

$$\mathrm{FC}\;\left(\mathrm{Fold}\;\mathrm{Change}\right)=\frac{\mathrm{norm}\;\left(\mathrm{THM}\right)}{\mathrm{norm}\;\left(\mathrm{WT}\right)}$$

Data standardization makes data comparable across different samples:

$$\mathrm{Norm}\;={\mathrm{reads}}_{\mathrm{num}}\;\times\frac{1000000}{{\mathrm{clean}}_{\mathrm{reads}}}$$

readsnum: Number of reads enriched in target region.

cleanreads: The total number of Reads obtained by filtering Raw Reads.

log2FC: The ratio of norm between two samples (groups) was taken as the logarithm base 2. An absolute value greater than 1 indicates a difference.

Data availability

All data generated or analysed during this study are included in this published article, its supplementary information files and publicly available repositories. The analysed Perl scripts and datasets are available in the Figshare repository, https://doi.org/https://doi.org/10.6084/m9.figshare.29143613 [82]. The genome data used in this study include three genomes generated in our previous work, which have been deposited in the Genome Warehouse of the National Genomics Data Center, China (https://ngdc.cncb.ac.cn/gwh), under accession numbers GWHBJYT00000000, GWHBJYR00000000, and GWHBJYS00000000 [83]. In addition, 35 genomes were obtained from the NCBI database (accession numbers are listed in Additional file 1: Table S1). The publicly available ChIP-seq and RNA-seq datasets of mouse limb buds at different developmental stages analysed in this study were obtained from the ENCODE database (https://www.encodeproject.org/). The corresponding accession numbers are listed in Additional file 1: Table S2.

Abbreviations

CNEs:

Conserved non-coding elements

TFs:

Transcription factors

GO:

Gene ontology

WT:

Wild-type mice

THM:

Transgenic homozygous mice

P1:

Postnatal day 1

DEGs:

Differentially expressed genes

KEGG:

Kyoto Encyclopedia of Genes and Genomes

FC:

Fold change

FDR:

False discovery rate

AER:

Apical ectodermal ridge

ZPA:

Zone of polarizing activity

GCPS:

Greig cephalopolysyndactyly syndrome

PHS:

Pallister-Killian syndrome

Hi-C:

High-throughput chromosome conformation capture

Evo-Devo:

Evolutionary biology and developmental biology

BF:

Bayes factor

PEI:

Polyethylenimine Linear

DMEM:

Dulbecco’s Modified Eagle Medium

RT:

Room temperature

PAM:

Protospacer adjacent motif

ICR:

Institute of Cancer Research mice

PCR:

Polymerase chain reaction

NJNU:

Nanjing Normal University

ARC:

Animal Resources Center

CUT&Tag:

Cleavage Under Targets and Tagmentation

RT-qPCR:

Reverse Transcription Quantitative Polymerase Chain Reaction

References

  1. Seki R, Kamiyama N, Tadokoro A, Nomura N, Tsuihiji T, Manabe M, Tamura K. Evolutionary and developmental aspects of avian-specific traits in limb skeletal pattern. Zoolog Sci. 2012;29(10):631–44.

    Article  PubMed  Google Scholar 

  2. Mann A, Pardo JD, Maddin HC. Snake-like limb loss in a Carboniferous amniote. Nat Ecol Evol. 2022;6(5):614–21.

    Article  PubMed  Google Scholar 

  3. Schwaner MJ, Lin DC, McGowan CP. Jumping mechanics of desert kangaroo rats. J Exp Biol. 2018;221(Pt 22):jeb186700.

    Article  PubMed  Google Scholar 

  4. Legreneur P, Monteil KM, Pelle E, Montuelle S, Bels V. Submaximal leaping in the grey mouse lemur. Zoology (Jena). 2011;114(4):247–54.

    Article  PubMed  Google Scholar 

  5. Lambert O, Bianucci G, Salas-Gismondi R, Di Celma C, Steurbaut E, Urbina M, de Muizon C. An amphibious whale from the Middle Eocene of Peru reveals early south pacific dispersal of quadrupedal cetaceans. Curr Biol. 2019;29(8):1352-1359 e1353.

    Article  CAS  PubMed  Google Scholar 

  6. Thewissen JG, Cooper LN, Clementz MT, Bajpai S, Tiwari BN. Whales originated from aquatic artiodactyls in the Eocene epoch of India. Nature. 2007;450(7173):1190–4.

    Article  CAS  PubMed  Google Scholar 

  7. Swank S, Sanger TJ, Stuart YE. (Non) parallel developmental mechanisms in vertebrate appendage reduction and loss. Ecol Evol. 2021;11(22):15484–97.

    Article  PubMed  PubMed Central  Google Scholar 

  8. Cooper LN, Sears KE, Armfield BA, Kala B, Hubler M, Thewissen JGM. Review and experimental evaluation of the embryonic development and evolutionary history of flipper development and hyperphalangy in dolphins (Cetacea: Mammalia). Genesis. 2018;56(1):e23076.

    Article  Google Scholar 

  9. Thewissen JG, Cohn MJ, Stevens LS, Bajpai S, Heyning J, Horton WE Jr. 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.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  10. Li J, Shang S, Fang N, Zhu Y, Zhang J, Irwin DM, Zhang S, Wang Z. Accelerated evolution of limb-related gene Hoxd11 in the common ancestor of cetaceans and ruminants (Cetruminantia). G3 (Bethesda). 2020;10(2):515–24.

    Article  CAS  PubMed  Google Scholar 

  11. Kim DW, Choi SH, Kim RN, Kim SH, Paik SG, Nam SH, Kim DW, Kim A, Kang A, Park HS. Comparative genomic analysis of the false killer whale (Pseudorca crassidens) LMBR1 locus. Genome. 2010;53(9):658–66.

    Article  CAS  PubMed  Google Scholar 

  12. Wang Z, Yuan L, Rossiter SJ, Zuo X, Ru B, Zhong H, Han N, Jones G, Jepson PD, Zhang S. Adaptive evolution of 5’HoxD genes in the origin and diversification of the cetacean flipper. Mol Biol Evol. 2009;26(3):613–22.

    Article  CAS  PubMed  Google Scholar 

  13. Sokolowski DJ, Vasquez OE, Wilson MD, Sokolowski MB, Anreiter I. Transcriptomic effects of the foraging gene shed light on pathways of pleiotropy and plasticity. Ann N Y Acad Sci. 2023;1526(1):99–113.

    Article  CAS  PubMed  Google Scholar 

  14. Gottschalk A, Sczakiel HL, Hulsemann W, Schwartzmann S, Abad-Perez AT, Grunhagen J, Ott CE, Spielmann M, Horn D, Mundlos S, et al. HOXD13-associated synpolydactyly: extending and validating the genotypic and phenotypic spectrum with 38 new and 49 published families. Genet Med. 2023;25(11):100928.

    Article  CAS  PubMed  Google Scholar 

  15. Xu F, Shangguan X, Pan J, Yue Z, Shen K, Ji Y, Zhang W, Zhu Y, Sha J, Wang Y, et al. HOXD13 suppresses prostate cancer metastasis and BMP4-induced epithelial-mesenchymal transition by inhibiting SMAD1. Int J Cancer. 2021;148(12):3060–70.

    Article  CAS  PubMed  Google Scholar 

  16. Roscito JG, Sameith K, Kirilenko BM, Hecker N, Winkler S, Dahl A, Rodrigues MT, Hiller M. Convergent and lineage-specific genomic differences in limb regulatory elements in limbless reptile lineages. Cell Rep. 2022;38(3):110280.

    Article  CAS  PubMed  Google Scholar 

  17. Chai S, Chong Y, Yin D, Qiu Q, Xu S, Yang G. Genomic insights into adaptation to bipedal saltation and desert-like habitats of jerboas. Sci China Life Sci. 2024;67(9):2003–15.

    Article  CAS  PubMed  Google Scholar 

  18. Roscito JG, Sameith K, Parra G, Langer BE, Petzold A, Moebius C, Bickle M, Rodrigues MT, Hiller M. Phenotype loss is associated with widespread divergence of the gene regulatory landscape in evolution. Nat Commun. 2018;9(1):4737.

    Article  PubMed  PubMed Central  Google Scholar 

  19. Hornblad A, Bastide S, Langenfeld K, Langa F, Spitz F. Dissection of the Fgf8 regulatory landscape by in vivo CRISPR-editing reveals extensive intra- and inter-enhancer redundancy. Nat Commun. 2021;12(1):439.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  20. Taskiran II, Spanier KI, Dickmanken H, Kempynck N, Pancikova A, Eksi EC, Hulselmans G, Ismail JN, Theunis K, Vandepoel R, et al. Cell-type-directed design of synthetic enhancers. Nature. 2024;626(7997):212–20.

    Article  CAS  PubMed  Google Scholar 

  21. Arnold M, Stengel KR. Emerging insights into enhancer biology and function. Transcription. 2023;14(1–2):68–87.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  22. Reilly SK, Yin J, Ayoub AE, Emera D, Leng J, Cotney J, Sarro R, Rakic P, Noonan JP. Evolutionary genomics. Evolutionary changes in promoter and enhancer activity during human corticogenesis. Science. 2015;347(6226):1155–9.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  23. Liu X, Zhang Y, Liu W, Li Y, Pan J, Pu Y, Han J, Orlando L, Ma Y, Jiang L. A single-nucleotide mutation within the TBX3 enhancer increased body size in Chinese horses. Curr Biol. 2022;32(2):480-487e486.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  24. Wang Z, Peng C, Wu W, Yan C, Lv Y, Li JT. Developmental regulation of conserved non-coding element evolution provides insights into limb loss in squamates. Sci China Life Sci. 2023;66(10):2399–414.

    Article  CAS  PubMed  Google Scholar 

  25. Kvon EZ, Kamneva OK, Melo US, Barozzi I, Osterwalder M, Mannion BJ, Tissieres V, Pickle CS, Plajzer-Frick I, Lee EA, et al. Progressive loss of function in a limb enhancer during snake evolution. Cell. 2016;167(3):633-642e611.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  26. Leal F, Cohn MJ. Loss and re-emergence of legs in snakes by modular evolution of sonic hedgehog and HOXD enhancers. Curr Biol. 2016;26(21):2966–73.

    Article  CAS  PubMed  Google Scholar 

  27. Strahle U, Rastegar S. Conserved non-coding sequences and transcriptional regulation. Brain Res Bull. 2008;75(2–4):225–30.

    Article  PubMed  Google Scholar 

  28. Paparidis Z, Abbasi AA, Malik S, Goode DK, Callaway H, Elgar G, deGraaff E, Lopez-Rios J, Zeller R, Grzeschik KH. Ultraconserved non-coding sequence element controls a subset of spatiotemporal GLI3 expression. Dev Growth Differ. 2007;49(6):543–53.

    Article  CAS  PubMed  Google Scholar 

  29. Creyghton MP, Cheng AW, Welstead GG, Kooistra T, Carey BW, Steine EJ, Hanna J, Lodato MA, Frampton GM, Sharp PA, et al. Histone H3K27ac separates active from poised enhancers and predicts developmental state. Proc Natl Acad Sci U S A. 2010;107(50):21931–6.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  30. Uebbing S, Kocher AA, Baumgartner M, Ji Y, Bai S, Xing X, Nottoli T, Noonan JP. Evolutionary innovations in conserved regulatory elements associate with developmental genes in mammals. Mol Biol Evol. 2024;41(10):msae199.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  31. Zhu C, Li S, Zhang D, Zhang J, Wang G, Zhou B, Zheng J, Xu W, Wang Z, Gao X, et al. Convergent degenerated regulatory elements associated with limb loss in limbless amphibians and reptiles. Mol Biol Evol. 2024;41(11):msae239.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  32. Cohn MJ, Tickle C. Developmental basis of limblessness and axial patterning in snakes. Nature. 1999;399(6735):474–9.

    Article  CAS  PubMed  Google Scholar 

  33. Rahmanl TMZ. Morphogenesis of the rudimentary hind-limb of the glass snake (Ophisaurus apodus Pallas). Development. 1974;32(2):431–43.

    Article  Google Scholar 

  34. Nemec S, Luxey M, Jain D, Huang Sung A, Pastinen T, Drouin J. Pitx1 directly modulates the core limb development program to implement hindlimb identity. Development. 2017;144(18):3325–35.

    Article  CAS  PubMed  Google Scholar 

  35. Rosin JM, Abassah-Oppong S, Cobb J. Comparative transgenic analysis of enhancers from the human SHOX and mouse Shox2 genomic regions. Hum Mol Genet. 2013;22(15):3063–76.

    Article  CAS  PubMed  Google Scholar 

  36. Wade C, Brinas I, Welfare M, Wicking C, Farlie PG. Twist2 contributes to termination of limb bud outgrowth and patterning through direct regulation of Grem1. Dev Biol. 2012;370(1):145–53.

    Article  CAS  PubMed  Google Scholar 

  37. Reilly SK, Noonan JP. Evolution of gene regulation in humans. Annu Rev Genomics Hum Genet. 2016;17:45–67.

    Article  CAS  PubMed  Google Scholar 

  38. Agrawal P, Heimbruch KE, Rao S. Genome-wide maps of transcription regulatory elements and transcription enhancers in development and disease. Compr Physiol. 2018;9(1):439–55.

    Article  PubMed  PubMed Central  Google Scholar 

  39. Malkmus J, Ramos Martins L, Jhanwar S, Kircher B, Palacio V, Sheth R, Leal F, Duchesne A, Lopez-Rios J, Peterson KA, et al. Spatial regulation by multiple Gremlin1 enhancers provides digit development with cis-regulatory robustness and evolutionary plasticity. Nat Commun. 2021;12(1):5557.

    Article  PubMed  PubMed Central  Google Scholar 

  40. Shen Y, Yue F, McCleary DF, Ye Z, Edsall L, Kuan S, Wagner U, Dixon J, Lee L, Lobanenkov VV, et al. A map of the cis-regulatory sequences in the mouse genome. Nature. 2012;488(7409):116–20.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  41. Woolfe A, Goodson M, Goode DK, Snell P, McEwen GK, Vavouri T, Smith SF, North P, Callaway H, Kelly K, et al. Highly conserved non-coding sequences are associated with vertebrate development. PLoS Biol. 2005;3(1):e7.

    Article  PubMed  Google Scholar 

  42. Liang N, Deme L, Kong Q, Sun L, Cao Y, Wu T, Huang X, Xu S, Yang G. Divergence of Tbx4 hindlimb enhancer HLEA underlies the hindlimb loss during cetacean evolution. Genomics. 2022;114(2):110292.

    Article  CAS  PubMed  Google Scholar 

  43. Infante CR, Park S, Mihala AG, Kingsley DM, Menke DB. Pitx1 broadly associates with limb enhancers and is enriched on hindlimb cis-regulatory elements. Dev Biol. 2013;374(1):234–44.

    Article  CAS  PubMed  Google Scholar 

  44. Ito S, Kitazawa R, Haraguchi R, Kondo T, Ouchi A, Ueda Y, Kitazawa S. Novel GLI3 variant causing overlapped Greig cephalopolysyndactyly syndrome (GCPS) and Pallister-Hall syndrome (PHS) phenotype with agenesis of gallbladder and pancreas. Diagn Pathol. 2018;13(1):1.

    Article  PubMed  PubMed Central  Google Scholar 

  45. Kuijper S, Feitsma H, Sheth R, Korving J, Reijnen M, Meijlink F. Function and regulation of Alx4 in limb development: complex genetic interactions with Gli3 and Shh. Dev Biol. 2005;285(2):533–44.

    Article  CAS  PubMed  Google Scholar 

  46. Osterwalder M, Barozzi I, Tissieres V, Fukuda-Yuzawa Y, Mannion BJ, Afzal SY, Lee EA, Zhu Y, Plajzer-Frick I, Pickle CS, et al. Enhancer redundancy provides phenotypic robustness in mammalian development. Nature. 2018;554(7691):239–43.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  47. Ali S, Abrar M, Hussain I, Batool F, Raza RZ, Khatoon H, Zoia M, Visel A, Shubin NH, Osterwalder M, et al. Identification of ancestral gnathostome Gli3 enhancers with activity in mammals. Dev Growth Differ. 2024;66(1):75–88.

    Article  CAS  PubMed  Google Scholar 

  48. Abbasi AA, Paparidis Z, Malik S, Bangs F, Schmidt A, Koch S, Lopez-Rios J, Grzeschik KH. Human intronic enhancers control distinct sub-domains of Gli3 expression during mouse CNS and limb development. BMC Dev Biol. 2010;10:44.

    Article  PubMed  PubMed Central  Google Scholar 

  49. Bolt CC, Lopez-Delisle L, Mascrez B, Duboule D. Mesomelic dysplasias associated with the HOXD locus are caused by regulatory reallocations. Nat Commun. 2021;12(1):5013.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  50. Marinic M, Aktas T, Ruf S, Spitz F. An integrated holo-enhancer unit defines tissue and gene specificity of the Fgf8 regulatory landscape. Dev Cell. 2013;24(5):530–42.

    Article  CAS  PubMed  Google Scholar 

  51. Kvon EZ, Waymack R, Gad M, Wunderlich Z. Enhancer redundancy in development and disease. Nat Rev Genet. 2021;22(5):324–36.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  52. Purushothaman S, Elewa A, Seifert AW. Fgf-signaling is compartmentalized within the mesenchyme and controls proliferation during salamander limb development. Elife. 2019;8:e48507.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  53. Li X, Cao X. BMP signaling and HOX transcription factors in limb development. Front Biosci. 2003;8:s805-812.

    Article  CAS  PubMed  Google Scholar 

  54. Geetha-Loganathan P, Nimmagadda S, Scaal M. Wnt signaling in limb organogenesis. Organogenesis. 2008;4(2):109–15.

    Article  PubMed  PubMed Central  Google Scholar 

  55. Zhulyn O, Li D, Deimling S, Vakili NA, Mo R, Puviindran V, Chen MH, Chuang PT, Hopyan S, Hui CC. A switch from low to high Shh activity regulates establishment of limb progenitors and signaling centers. Dev Cell. 2014;29(2):241–9.

    Article  CAS  PubMed  Google Scholar 

  56. 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.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  57. Chen C, Chen H, Zhang Y, Thomas HR, Frank MH, He Y, Xia R. TBtools: an integrative toolkit developed for interactive analyses of big biological data. Mol Plant. 2020;13(8):1194–202.

    Article  CAS  PubMed  Google Scholar 

  58. Kumar S, Stecher G, Suleski M, Hedges SB. TimeTree: a resource for timelines, timetrees, and divergence times. Mol Biol Evol. 2017;34(7):1812–9.

    Article  CAS  PubMed  Google Scholar 

  59. Nguyen LT, Schmidt HA, von Haeseler A, Minh BQ. IQ-TREE: a fast and effective stochastic algorithm for estimating maximum-likelihood phylogenies. Mol Biol Evol. 2015;32(1):268–74.

    Article  CAS  PubMed  Google Scholar 

  60. Hecker N, Hiller M. A genome alignment of 120 mammals highlights ultraconserved element variability and placenta-associated enhancers. Gigascience. 2020;9(1):giz159.

    Article  PubMed  PubMed Central  Google Scholar 

  61. Blanchette M, Kent WJ, Riemer C, Elnitski L, Smit AF, Roskin KM, Baertsch R, Rosenbloom K, Clawson H, Green ED, et al. Aligning multiple genomic sequences with the threaded blockset aligner. Genome Res. 2004;14(4):708–15.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  62. Hubisz MJ, Pollard KS, Siepel A. PHAST and RPHAST: phylogenetic analysis with space/time models. Brief Bioinform. 2011;12(1):41–51.

    Article  CAS  PubMed  Google Scholar 

  63. Hu Z, Sackton TB, Edwards SV, Liu JS. Bayesian detection of convergent rate changes of conserved noncoding elements on phylogenetic trees. Mol Biol Evol. 2019;36(5):1086–100.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  64. McLean CY, Bristor D, Hiller M, Clarke SL, Schaar BT, Lowe CB, Wenger AM, Bejerano G. GREAT improves functional interpretation of cis-regulatory regions. Nat Biotechnol. 2010;28(5):495–501.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  65. Shen WK, Chen SY, Gan ZQ, Zhang YZ, Yue T, Chen MM, Xue Y, Hu H, Guo AY. AnimalTFDB 4.0: a comprehensive animal transcription factor database updated with variation and expression annotations. Nucleic Acids Res. 2023;51(D1):D39–45.

    Article  CAS  PubMed  Google Scholar 

  66. Lambert SA, Yang AWH, Sasse A, Cowley G, Albu M, Caddick MX, Morris QD, Weirauch MT, Hughes TR. Similarity regression predicts evolution of transcription factor sequence specificity. Nat Genet. 2019;51(6):981–9.

    Article  CAS  PubMed  Google Scholar 

  67. Zhu B, Yang C, Sun L, Li Z, Li J, Hua ZC. Expression pattern and prognostic implication of zinc homeostasis-related genes in acute myeloid leukemia. Metallomics. 2023;15(5):mfad022.

    Article  PubMed  Google Scholar 

  68. Rauluseviciute I, Riudavets-Puig R, Blanc-Mathieu R, Castro-Mondragon JA, Ferenc K, Kumar V, Lemma RB, Lucas J, Cheneby J, Baranasic D, et al. JASPAR 2024: 20th anniversary of the open-access database of transcription factor binding profiles. Nucleic Acids Res. 2024;52(D1):D174–82.

    Article  CAS  PubMed  Google Scholar 

  69. Bhattacharya P, Edwards K, Schmid KL. Segmentation methods and morphometry of confocal microscopy imaged corneal epithelial cells. Cont Lens Anterior Eye. 2022;45(6):101720.

    Article  PubMed  Google Scholar 

  70. Booth M, Powell N, Corfield C, French JM. An automated technique for double staining of bone and cartilage in fetal mouse skeletal specimens using alizarin red S and Alcian blue. Biotech Histochem. 2022;97(3):222–7.

    Article  CAS  PubMed  Google Scholar 

  71. Patuel SJ, English C, Lopez-Scarim V, Konig I, Souders CL 2nd, Ivantsova E, Martyniuk CJ. Dataset for diseases associated with exposure to broflanilide, a novel pesticide, in larval zebrafish (Danio rerio). Data Brief. 2023;50:109534.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  72. Anders S, Huber W. Differential expression analysis for sequence count data. Genome Biol. 2010;11(10):R106.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  73. Conesa A, Gotz S, Garcia-Gomez JM, Terol J, Talon M, Robles M. Blast2GO: a universal tool for annotation, visualization and analysis in functional genomics research. Bioinformatics. 2005;21(18):3674–6.

    Article  CAS  PubMed  Google Scholar 

  74. Tatusov RL, Galperin MY, Natale DA, Koonin EV. The COG database: a tool for genome-scale analysis of protein functions and evolution. Nucleic Acids Res. 2000;28(1):33–6.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  75. Kang Y, Kim YW, Kang J, Kim A. Histone H3K4me1 and H3K27ac play roles in nucleosome eviction and eRNA transcription, respectively, at enhancers. FASEB J. 2021;35(8):e21781.

    Article  CAS  PubMed  Google Scholar 

  76. Kaya-Okur HS, Wu SJ, Codomo CA, Pledger ES, Bryson TD, Henikoff JG, Ahmad K, Henikoff S. CUT&Tag for efficient epigenomic profiling of small samples and single cells. Nat Commun. 2019;10(1):1930.

    Article  PubMed  PubMed Central  Google Scholar 

  77. Bolger AM, Lohse M, Usadel B. Trimmomatic: a flexible trimmer for Illumina sequence data. Bioinformatics. 2014;30(15):2114–20.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  78. Li H, Durbin R. Fast and accurate short read alignment with Burrows-Wheeler transform. Bioinformatics. 2009;25(14):1754–60.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  79. Li H, Handsaker B, Wysoker A, Fennell T, Ruan J, Homer N, Marth G, Abecasis G, Durbin R, Genome Project Data Processing S. The sequence alignment/map format and SAMtools. Bioinformatics. 2009;25(16):2078–9.

    Article  PubMed  PubMed Central  Google Scholar 

  80. Grytten I, Rand KD, Nederbragt AJ, Storvik GO, Glad IK, Sandve GK. Graph peak caller: calling ChIP-seq peaks on graph-based reference genomes. PLoS Comput Biol. 2019;15(2):e1006731.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  81. Love MI, Huber W, Anders S. Moderated estimation of fold change and dispersion for RNA-seq data with DESeq2. Genome Biol. 2014;15(12):550.

    Article  PubMed  PubMed Central  Google Scholar 

  82. Zhang Z, Zhenpeng Yu, Chong Y, Liu Y, Liu J, Ren W, Shixia Xu, Yang G. Evolution of cetacean-specific conserved non-coding elements suggests their role in the limb changes during secondary aquatic adaptation. 2025. FigShare. https://doi.org/10.6084/m9.figshare.29143613.

  83. Xu S, Shan L, Tian R, Yu Z, Sun D, Zhang Z, Seim I, Zhou M, Sun L, Liang N, et al. Multi-level genomic convergence of secondary aquatic adaptation in marine mammals. Innovation (Camb). 2025;6(3):100798.

    CAS  PubMed  Google Scholar 

Download references

Acknowledgements

We thank all the members of Jiangsu Key Laboratory for the Conservation and Utilization of Biodiversity in the Middle and Lower Reaches of the Yangtze River for the daily discussion.

Funding

This research was financially supported by the Key Project of the National Natural Science Foundation of China (grant no. 32030011), the National Key Research and Development Program, Ministry of Science and Technology (grant no. 2022YFF1301600), and the National Natural Science Foundation of China (grant nos. 32070409, 32270453), the Priority Academic Program Development of Jiangsu Higher Education Institutions (PAPD), the Qinglan Project of Jiangsu Province, and the PI Project of Southern Marine Science and Engineering Guangdong Laboratory (Guangzhou) (grant no. GML2021GD0805).

Author information

Authors and Affiliations

  1. Jiangsu Key Laboratory for the Conservation and Utilization of Biodiversity in the Middle and Lower Reaches of the Yangtze River, College of Life Sciences, Nanjing Normal University, Nanjing, 210023, China

    Zhenhua Zhang, Zhenpeng Yu, Yujie Chong, Yao Liu, Jia Liu, Wenhua Ren, Shixia Xu & Guang Yang

  2. Southern Marine Science and Engineering Guangdong Laboratory (Guangzhou), Guangzhou, 511458, China

    Guang Yang

Authors
  1. Zhenhua Zhang
    View author publications

    Search author on:PubMed Google Scholar

  2. Zhenpeng Yu
    View author publications

    Search author on:PubMed Google Scholar

  3. Yujie Chong
    View author publications

    Search author on:PubMed Google Scholar

  4. Yao Liu
    View author publications

    Search author on:PubMed Google Scholar

  5. Jia Liu
    View author publications

    Search author on:PubMed Google Scholar

  6. Wenhua Ren
    View author publications

    Search author on:PubMed Google Scholar

  7. Shixia Xu
    View author publications

    Search author on:PubMed Google Scholar

  8. Guang Yang
    View author publications

    Search author on:PubMed Google Scholar

Contributions

G.Y, W.R and S.X conceptualized and supervised the study. Z.Z performed all sample collection, phenotype and data analysis, RNA extraction, and RT-PCR validation. Z.Z and Z.Y were responsible for the bioinformatic analysis. Z.Z, Y.L, J.L, and Y.C were responsible for mouse husbandry and genotype identification. Z.Z prepared the original draft, and G.Y revised the manuscript. All authors read and approved the final manuscript.

Corresponding authors

Correspondence to Shixia Xu or Guang Yang.

Ethics declarations

Ethics approval and consent to participate

All animal experiments in this manuscript were conducted according to the Animal Care and Use Committee of Nanjing Normal University (ethics approval number: IACUC-20200501).

Consent for publication

Not applicable.

Competing interests

The authors declare no competing interests.

Additional information

Publisher’s Note

Springer Nature remains neutral with regard to jurisdictional claims in published maps and institutional affiliations.

Supplementary Information

12915_2025_2300_MOESM1_ESM.xlsx

Additional file 1. Table S1. List of species used to identify specific changes in CNEs in cetaceans. Table S2 Public data from the ENCODE database. Table S3 Sequence-divergent CNEs in cetaceans associated with limb development. Table S4 Functional annotation of cetacean-specific CNEs was performed using GREAT. Table S5 163 cetacean-specific CNEs potentially associated with cetacean limb changes. Table S6 AnimalTFDB-4.0 predicts transcription factors that might bind mouse CNEs and cetacean CNEs, respectively. Table S7 CIS-BP predicts transcription factors that might bind mouse CNEs and cetacean CNEs, respectively. Table S8 To analyze the co-expression between the predicted transcription factors and target genes based on Spearman. Table S9 Biological process enrichment analysis of transcription factors significantly associated with target gene expression was performed using Metascape. Table S10 Transcription factor predictions for cetacean-specific CNE InDel sites were performed based on the JASPAR public database. Table S11 Forelimb and Hindlimb bud area of WT and THM at E10.5–E11.5. Table S12 The normality of the morphological data of E10.5–E11.5 limb buds was tested based on skewness-kurtosis and the Shapiro-Wilk normality test. Table S13 Characterization of forelimb and hindlimb elements in WT and THM at E18.5. Table S14 The normality of the morphological data of E18.5 limb was tested based on Skewness-Kurtosis and Shapiro-Wilk normality test. Table S15 Characterization of limb elements in wild-type and hs1586 knock-in mice on postnatal one day. Table S16 The normality of the morphological data of postnatal one-day limbs was tested based on Skewness-Kurtosis and Shapiro-Wilk normality test. Table S17 Differentially expressed genes (DEGs) of the hs1586 knock in (KI) mice compared with the wide-type (WT) embryos at E10.5. Table S18 Differentially expressed genes (DEGs) of the hs1586 knock in (KI) mice compared with the wide-type (WT) embryos at E11.5. Table S19 Differentially expressed genes (DEGs) of the hs1586 knock in (KI) mice compared with the wide-type (WT) embryos at E12.5. Table S20 CUT&Tag analysis of the cetaceans hs1586 knock in (KI) mice and wild type embryos hindlimb bub at E10.5-E11.5. Table S21 Functional analysis validated conserved non-coding elements and enhancers around the mouse Gli3 gene (mm10). Table S22 Primer sequences for real-time quantitative RT-PCR and hs1586 sequence of cetacean and mouse in this study.

12915_2025_2300_MOESM2_ESM.pdf

Additional file 2. Fig. S1 Collinearity analysis of the genomes and filtration of CNEs. Fig. S2 Functional enrichment analysis of cetacean sequence-divergent CNEs associated with limb development in GO biological processes was performed using GREAT. Fig. S3 Functional enrichment analysis of cetacean sequence-divergent CNEs associated with limb development in mammalian phenotypes was performed using GREAT. Fig. S4 Filter CNEs and select cetacean-specific CNEs closely associated with limb development. Fig. S5 The nucleotide substitution rate of DNA was detected using PhyloAcc software (see Materials and Methods), CNE90, CNE227, CNE187, and CNE239 exhibited the signal of accelerated evolution in cetaceans, and these CNEs are associated with genes of Pitx1, Bmp2, Tcf4, and Tbx15, respectively. Blue represents cetaceans (baleen whales and toothed whales). Fig. S6 Using the PhyloAcc software (see Methods), we detected signals of accelerated evolution in cetaceans for CNE320 and CNE315 associated with Tbx3, as well as CNE218 and CNE219 associated with Alx4. Blue represents cetaceans (baleen whales and toothed whales). Fig. S7 Based on multi-genome sequence alignment, CNE622, CNE682, CNE497, CNE531, and CNE629 exhibit species-specific deletions in cetaceans (Blue dashed box), and these CNEs are associated with genes of Shox2, Hoxa13, Pax9, Wnt5a, and Tbx15, respectively. Fig. S8 The Gli3 enhancers hs1586 and CNE102 (mm1887) exhibit sequence divergence in cetaceans. Fig. S9 Relative luciferase activity (Fluc/Rluc) of eleven cetacean-specific CNEs in Control, M. musculus, B. acutorostrata, and T. truncatus groups. a-k. Data are mean±standard error of three individual experiments. P-value was calculated using Student’s t-test (*, p<0.05; **, p<0.01; ***, p<0.001). Fig. S10 Generation and identification of a transgenic mouse model carrying the cetacean-specific enhancer hs1586. Fig. S11 Morphometric analyses of the limb buds in WT and THM at E10.5. Fig. S12 Morphometric analyses of the limb buds in WT and THM at E11.5. Fig. S13 Morphometric analyses of the limb buds in WT and THM at E18.5 and postnatal day one. Fig. S14 Before RNA-seq analysis, THM were confirmed using PCR and agarose gel electrophoresis. M: blank control, C: positive control (wild-type), W: wild-type, T: transgenic homozygous mice. Fig. S15 RNA-seq analysis of hindlimb at E10.5. Fig. S16 RNA-seq analysis of Gli3 and Inhba

Rights and permissions

Open Access This article is licensed under a Creative Commons Attribution-NonCommercial-NoDerivatives 4.0 International License, which permits any non-commercial use, sharing, distribution and reproduction in any medium or format, as long as you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons licence, and indicate if you modified the licensed material. You do not have permission under this licence to share adapted material derived from this article or parts of it. The images or other third party material in this article are included in the article’s Creative Commons licence, unless indicated otherwise in a credit line to the material. If material is not included in the article’s Creative Commons licence and your intended use is not permitted by statutory regulation or exceeds the permitted use, you will need to obtain permission directly from the copyright holder. To view a copy of this licence, visit http://creativecommons.org/licenses/by-nc-nd/4.0/.

Reprints and permissions

About this article

Check for updates. Verify currency and authenticity via CrossMark

Cite this article

Zhang, Z., Yu, Z., Chong, Y. et al. Evolution of cetacean-specific conserved non-coding elements suggests their role in the limb changes during secondary aquatic adaptation. BMC Biol 23, 216 (2025). https://doi.org/10.1186/s12915-025-02300-0

Download citation

  • Received:

  • Accepted:

  • Published:

  • Version of record:

  • DOI: https://doi.org/10.1186/s12915-025-02300-0

Keywords

BMC Biology

ISSN: 1741-7007

Contact us