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

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.

Keywords: CNEs; Cetaceans; Limb morphology; Phenotypic recovery; Transgenic mice.

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

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

Fig. 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)

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

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