Developmental hourglass and heterochronic shifts in fin and limb development

  1. Koh Onimaru  Is a corresponding author
  2. Kaori Tatsumi
  3. Chiharu Tanegashima
  4. Mitsutaka Kadota
  5. Osamu Nishimura
  6. Shigehiro Kuraku  Is a corresponding author
  1. Laboratory for Phyloinformatics, RIKEN Center for Biosystems Dynamics Research (BDR), Japan
  2. Laboratory for Bioinformatics Research, RIKEN BDR, Japan
  3. Molecular Oncology Laboratory, Graduate School of Medicine, Nagoya University, Japan
6 figures, 3 tables and 12 additional files


Figure 1 with 8 supplements
Transcriptome analysis and orthology assignment.

(A) The skeletal patterns of a mouse limb (top) and a bamboo shark pectoral fin (bottom). Anterior is to the top; distal is to the right. (B) Mouse forelimb buds and bamboo shark pectoral fin buds that were analyzed by RNA-seq. (C) Comparison of the accuracy of three orthology assignment methods. Vertical axis, the percentages of correctly assigned Hoxa and Hoxd paralogs (black bars) and Fgf paralogs (white bars). (D) Heat map visualization of the transcription profile of Hoxa and Hoxd genes in mouse limb buds (left) and bamboo shark fin buds (right) with scaled TPMs.

Figure 1—figure supplement 1
Schematic representation of the orthology assignment algorithm.

Red arrows, the main flow of the algorithm. Black arrows, orthology assignment for cartilaginous fish-specific genes. Gray arrows, parallel retrieving of orthologs of mouse genes from other animals. Red rectangles, best hits across other animal genes or in elephantfish genes. Green rectangles, best hits among each animal genome. Note that this schematic explains how the orthology of abstract genes ‘bamboo shark gene X’ and ‘mouse gene Y’ are assigned. First, using BLASTP, putative orthologs of bamboo shark genes are retrieved from other animal genomes, such as human, mouse, alligator, and elephantfish. Then, all BLASTP results except those from elephantfish are concatenated to find a best scored gene across species (cross-species best hit). In this schematic, the alligator gene XP001 is the best hit. In parallel, putative orthologs of mouse genes are also retrieved from the same set of animal genomes. If there is a mouse gene Y that has a best hit with alligator XP001, this mouse gene Y and bamboo shark gene X are considered to be an orthologous pair.

Figure 1—figure supplement 2
Molecular phylogenetic tree for Fgf family.

The tree was inferred with the maximum-likelihood method. The support values at nodes indicate bootstrap probabilities. Genes highlighted in red are bamboo shark genes (can be converted into the original gene ID by replacing ‘g’ with ‘Chipu’ and fill the digit with 0 to be seven figure number in total).

Figure 1—figure supplement 3
Additional molecular phylogenetic trees for Fgf8, Fgf11, Fgf12, and Fgf13.

These trees are shown because alignment sequences used in Figure 1—figure supplement 2 are truncated or absent in these genes. The tree was inferred with the maximum-likelihood method. The support values at nodes indicate bootstrap probabilities. Genes highlighted in red are bamboo shark genes.

Figure 1—figure supplement 4
Comparison between the TPM and TMM.

(A) Visualization of the effect of normalization by showing a housekeeping gene family, Ndufa. Left panels show TMM (trimmed mean of M) and TPM (transcripts per million) calculated by RSEM. Right panels show these values with additional normalization using several other housekeeping genes (Atp5j, Atp5h, Atp5g3, Psmc3, Psmc5, Psmd7, Mrpl54, Mrpl46, Polr2e, Polr1b, Mrpl2). Housekeeping genes are selected from a previously published list (; Eisenberg and Levanon, 2013). All expression values are standardized by setting the maximum expression value of each gene as ‘1’. Note that because housekeeping genes do not change their expression amount over time, these expression values should be close to ‘1’ (i.e. all colors should be dark blue) with some exceptions. However, the intact TMM (top, left panel) is apparently biased, in that the majority of Ndufa genes show their strongest expression at E9.5, with sharp decreases at other stages. This bias can be corrected by normalization with other housekeeping genes (top right panel). In contrast, the intact TPM (bottom, left panel) has a weaker bias than TMM. Additional normalization (bottom, right panel) has less of an effect. Therefore, this study used the intact TPM. (B) Euclidean distances of transcriptome data between mouse samples (left) and between bamboo shark samples (right). Whereas the close relation of the replicates of mouse samples can be seen from this heat map, the replicates of bamboo shark samples show less similarity. This noisy data may be attributed to the fact that there is no established strain of the bamboo shark and/or that bamboo shark embryos were staged by morphology but not physical time. However, the average of replicates seems to mitigate the noise of the bamboo shark samples, because Hox gene expression showed a smooth temporal collinearity as seen in Figure 1D.

Figure 1—figure supplement 5
The effect of scaling methods to housekeeping genes.

(A) A simple example for comparing expression distances between two species. Species 1 and 2 are imaginery simple species that have two genes (genes 1 and 2) and three developmental time points (t1, t2, and t3). Distances in the bottom are the Euclidean distance between two species at each stage. (B, C) Housekeeping gene expressions with intact TPM values and different scaling methods in mouse limb buds (B) and bambooshark fin buds (C). Intact TPM, TPM values without any scaling methods; Max 1, TPM values were scaled by setting the highest TPM in each gene of each species to ‘1’; Z-score, the mean expression value was subtracted from each expression value and each result was then divided by the standard deviation; Unit vector, expression values were divided by the norm; Log10, log10 transformation. These housekeeping genes are listed in both a human housekeeping gene list ( and the BUSCO data set (thus these genes are likely conserved in most vertebrates). Note that although the expression values of the housekeeping genes were almost constant during development, Z-score scaling amplifies subtle differences between stages. In addition, intact TPM values were not readily comparable between limb buds and bamboo shark fin buds (e.g. the maximum value of POLR1B in mouse limb buds was roughly twice as high as that of bamboo shark fin buds). Error bars are not displayed.

Figure 1—figure supplement 6
The effect of scaling methods to heterochronic genes.

(A, B) Heterochronic gene expressions with intact TPM values and different scaling methods in mouse limb buds (A) and bambooshark fin buds (B). See the legend of Figure 1—figure supplement 5 for scaling methods. Error bars are omitted. (C) The total Euclidean distance with respect to gene expression for the three housekeeping and the three heterochronic genes between mouse limb buds and bamboo shark fin buds. Using the housekeeping genes shown in (A, B) and the different scaling methods, the graph shows the summation of Euclidean distances between all combinations of mouse limb and bamboo shark fin stages. (D) The ratios of the Euclidean distances for housekeeping genes to those for heterochronic genes as shown in C.

Figure 1—figure supplement 7
Examination of quantitative collinearity of Hoxd genes.

TPM values of 5ʹ Hoxd genes in mouse limb buds at E12.5 (left) and bamboo shark fin buds at st. 31 (right, orange bars) and st. 32 (right, blue bars). Error bars, SEM. Note that the genomic locus of Hoxd genes was positively correlated with their expression amount in the mouse limb bud, whereas no such correlation was found in the bamboo shark fin bud.

Figure 1—figure supplement 8
Expression profile of genes related to cellular differentiation.

(A, B) Scaled TPM values of indicated genes related to chondrogenesis (A) and myogenesis (B). Error bars, SEM.

Figure 2 with 1 supplement
Detection of heterochronic gene expression between mouse limb buds and bamboo shark fin buds.

(A) Clustering analysis of gene expression dynamics. Each column represents an ortholog pair between the bamboo shark and the mouse. Each row indicates scaled gene expression at a time point indicated to the right of the heat map. Values are scaled TPMs. (B, C) Whole-mount in situ hybridization of Hand2 (B) and Vcan (C) as examples of the heterochronic genes detected in (A). Asterisks, background signals; scale bars, 200 μm. Error bars: SEM.

Figure 2—figure supplement 1
Other heterochronic genes.

(A) Left panels, whole mount in situ hybridization of Aldh1a2 (one of the genes from the cluster Heterochronic1) in bamboo shark fin buds and mouse limb buds. Right panels, TPM values of Aldh1a2. Arrowheads indicate the late-stage expression of Aldh1a2 in bamboo shark fin buds. Scale bars, 200 μm. Error bars, SEM. (B) Heatmap of genes that exhibit an inverse relation to Heterochronic2 genes in Figure 2A. Yellow empty box, genes that exhibit relatively sharp upregulation in bamboo shark fin buds and downregulation in mouse limb buds over time. (C) Comparison of Fgf gene expression. Vertical axis, TPM values; error bars, SEM; N/A, not applicable because of the absence of Fgf24 in the mouse genome.

Figure 3 with 1 supplement
Shh pathway in mouse limb buds and bamboo shark fin buds.

(A, B) Scaled expression of Shh and related genes in mouse limb buds (A) and bamboo shark fin buds (B), respectively. The rectangles indicate the expression peaks of Shh, Hoxd9, and Hoxd10 (magenta), Shh target genes (yellow) and Hoxd11 and Hoxd12 (green). (C, D) Whole-mount in situ hybridization of Ptch1 and Hoxd12 in mouse limb buds (C) and bamboo shark fin buds (D); scale bars, 200 μm. White arrowheads in (D) indicate restricted expression of Ptch1 in bamboo shark fin buds. Black arrowheads in (C and D) indicate anteriorly extended expression of Hoxd12.

Figure 3—figure supplement 1
The temporal dynamics of the Shh pathway based on intact TPM values.

Red filled rectangles, the expression peak of Shh, Hoxd9, and Hoxd10; yellow filled rectangles, the expression peak of Shh target genes; green filled rectangles, the expression peak of Hoxd11 and Hoxd12.

Figure 4 with 2 supplements
Hourglass-shaped conservation of the transcriptome profile between fins and limbs.

(A) Euclidean distances of the transcriptome profiles. Every combination of time points of bamboo shark fin buds and mouse limb buds is shown. The darker colors indicate a greater similarity between gene expression profiles. (B) A line plot of the Euclidean distances shown in (A). The x axis indicates the mouse limb stages, and the y axis is the Euclidean distance. (C) The same as (A) except that only Hoxd genes are included. (D, E) Scatter plots of the first and second principal components (D) and of the second and third components (e). Arrows in (E) indicate the time-order of transcriptome data. (F) Count of tissue-associated genes expressed in mouse forelimb buds. Genes with 0.65 ≤ entropy were counted.

Figure 4—figure supplement 1
Confirmation analyses of the transcriptome comparison.

Cross-species comparisons of transcriptome data between the two species with indicated distance methods. Note that these methods consistently show the closest distance around E10.5 and st. 27.5–30.

Figure 4—figure supplement 2
Additional PCA data and counts for stage- and tissue-associated genes.

(A) The ratio of explained variable for each of the principal components from Figure 4D and E. (B) Euclidean distance measures using the indicated principal components. Note that individual principal components do not reproduce the hourglass-shaped conservation shown in Figure 4A, but PC1, PC2, and PC3 are sufficient for the most part to reproduce Figure 4A. (C) Percentage of stage-associated genes with [z ≤ 1.0] for mouse limb buds (left) and bamboo shark fin buds (right). Note that both species showed a low percentage of stage-associated genes during the middle stages of development. (D) Number (left) and fraction (right) of tissue-associated genes expressed in mouse limb buds. Tissue specificity was evaluated by entropy using RNA-seq data from 71 mouse tissues. A gene with entropy ≥0.65 was considered a tissue-specific gene. In the right panel, gene counts were normalized based on the number of total expressed genes. Note that the number of tissue-associated genes was lowest at E10.5.

Figure 5 with 2 supplements
Hourglass-shaped conservation of OCRs in mouse limb development.

(A) ATAC-seq signals in the enhancer regions of the HoxA cluster. e1 to e4, known limb enhancers. Green vertical lines below the signals, peak regions determined by MACS2. (B) Comparison of a quality index, FRiP, for ATAC-seq data. Blue bars are samples with a FRiP score >0.2. The number in the end of the label name indicates the replicate number. (C) Conservation analysis of sequences in ATAC-seq peaks with BLASTN. The values to the right of each graph indicate the fraction of conserved sequences in the total peak regions. The common name of each genome sequence is indicated above the graph. The not-conserved heatmap indicates the fraction of sequences that were not aligned to any genome sequences and thus serve as a negative control. (D) Temporal changes of sequence conservation frequency in ATAC-seq peaks with LAST. Error bars: SEM.

Figure 5—figure supplement 1
ATAC-seq quality control.

(A) Correlation distance between samples. The numbers in the end of the sample names indicate the replicates of indicated stages. Darker color means more similar gene expressions. (B) Percentage of peak regions in the genome sequence. (C) ATAC-seq signals in BPM (blue signals), peak regions (blue rectangles) and the known limb enhancers of HoxA cluster (red rectangles, e1–e19). Note that only e5 is not covered by ATAC-seq data.

Figure 5—figure supplement 2
Conservation measures of OCRs.

(A, B) The absolute count of OCRs that overlap with sequences conserved between the mouse and the alligator (A) and the bamboo shark (B). Error bars, SEM. (C, D) The fraction of conserved OCRs sorted by the identified clusters in Figure 6A. Sequence conservation was estimated by pairwise alignment using LAST (A–D).

Figure 6 with 5 supplements
Temporal dynamics of OCRs during mouse limb development.

(A) The heatmap (left) shows whole-genome clustering of ATAC-seq peaks. Each row indicates a particular genome region with a length of 1400 bp. Columns indicate developmental stages. C1−C8 are cluster numbers. The motifs (right) show the rank of enriched motifs in the sequences of each cluster. (B) Top, volcano plots of ATAC-seq signals between indicated stages (p-values, two-sided Student’s t-test). Bottom, the counts of differential signals (black dots in the top panel). + and − are genomic regions with increased or decreased signals, respectively. (C) The fraction of limb-specific OCRs for each cluster.

Figure 6—figure supplement 1
Clustering analyses of ATAC-seq peaks with different replicates.

DDifferent replicates were used for the same analysis as shown in Figure 6A. The number after the stage name indicates the replication number. Note that the clustering analyses with different replicates identified clusters similar to those in Figure 6A (compare the left-most panel with the second and third panels from the left). Including replicates with a low-quality score resulted in a relatively small fraction of early stage-specific peaks and a large fraction of late stage-specific peaks (right-most panel).

Figure 6—figure supplement 2
Analysis of enrichment for known motifs in ATAC-seq peaks.

The top five motifs from clusters C5 and C6 determined while using all other peaks as the background sequence.

Figure 6—figure supplement 3
De novo motif discoveries and known motif enrichment analysis of ATAC-seq peaks with an alternative background.

The top five motifs from clusters C5 and C6 determined while using all other peaks as the background sequence.

Figure 6—figure supplement 4
De novo motif discoveries of ATAC-seq peaks.

The top five motifs from each cluster. See Supplementary data for the full list of motifs. C1C8 correspond to the clusters in Figure 6A.

Figure 6—figure supplement 5
Counts of accessible motifs at each stage.

The average number of top-ranked motifs identified by de novo motif discovery in Figure 6—figure supplement 4 are plotted against mouse embryonic stages. The average numbers were calculated using all three replicates of ATAC-seq peaks at each stage. Rows indicate clusters identified in Figure 6A; columns indicate motif rank. Error bars, SEM. C1C8 correspond to the clusters in Figure 6A. Note that the number of CTCF motifs (top-ranked in C5) was relatively stable over time, which is consistent with the clustering analysis shown in Figure 6A. In addition, the number of motifs enriched for C3, such as BHLHA15, HOX13, TEAD, and Tlx?, increased over time. In contrast, COUP-TFII and TCF7L2 motifs decreased over time. Interestingly, LHX and HOX9 motifs were transiently increased at E10.5.


Table 1
Assembly statistics of bamboo shark transcriptome.
CharacteristicBamboo shark transcriptomeBamboo shark gene model (Hara et al., 2018)
Total number of sequences22201534038
Total sequence length (bp)19554136736633751
Average length (bp)8801076
Maximum length (bp)18451108594
N count010208
N50 length (bp)20751749
Protein coding6389834038
Orthology detected4163318180
Unique orthologs1413914907
Unique orthologs without gene symbols18211780
Unique orthologs only in elephantfish826552
Sequences with no orthology2089215254
Orthologs with mouse genes1232613005
Table 2
PCA loadings.
Loading axis: PC2
Gene symbolCluster nameLoading
Key resources table
Reagent type (species) or resourceDesignationSource or referenceIdentifiersAdditional information
Gene (Mus musculus)Hand2ENSEMBLENSMUST00000040104.4N/A
Gene (Mus musculus)VcanENSEMBLENSMUST00000109546.8N/A
Gene (Mus musculus)Aldh1a2ENSEMBLENSMUST00000034723.5N/A
Gene (Mus musculus)Ptch1ENSEMBLENSMUST00000192155.5N/A
Gene (Mus musculus)Hoxd12ENSEMBLENSMUST00000001878.5N/A
Gene (Chiloscyllium punctatum)Hand2ENSEMBLChipun0004250/g4250.t1/ TRINITY_DN85524_c0_g1_i1N/A
Gene (Chiloscyllium punctatum)VcanThis paperChipun0003941/g3941.t1/ TRINITY_DN95522_c0_g1_i8N/A
Gene (Chiloscyllium punctatum)Hoxd12This paperChipun0005654/g5654.t1/TRINITY_DN85970_c0_TRINITYg1_i1N/A
Gene (Chiloscyllium punctatum)Ptch1This paperChipun0003320/g3320.t1/TRINITY_DN92499_c0_g1_i3N/A
Gene (Chiloscyllium punctatum)Aldh1a2This paperChipun0010503/g10503.t1/TRINITY_DN81423_c0_g1_i1N/A
Strain, strain background (Mus musculus)C52BL/6Laboratory for Animal Resources and Genetic Engineering RIKEN,N/AN/A
AntibodyAnti-Digoxigenin-AP, Fab fragments (Sheep)Millipore SigmaCat# 11093274910polyclonal
Sequence-based reagentMus musculus Hand2 forward primerThis paperPCR primersACCAAACTCTCCAAGATCAAGACACTG
Sequence-based reagentMus musculus Hand2 reverse primerThis paperPCR primersTTGAATACTTACAATGTTTACACCTTC
Sequence-based reagentMus musculus Vcan forward primerThis paperPCR primersTGCAAAGATGGTTTCATTCAGCGACAC
Sequence-based reagentMus musculus Vcan reverse primerThis paperPCR primersACACGTGCAGAGACCTGCAAGATGCTG
Sequence-based reagentMus musculus Aldh1a2 forward primerThis paperPCR primersACCGTGTTCTCCAACGTCACTGATGAC
Sequence-based reagentMus musculus Aldh1a2 reverse primerThis paperPCR primersTCTGTCAGTAACAGTATGGAGAGCTTG
Sequence-based reagentMus musculus Hoxd12 forward primerThis paperPCR primersCTCAACTTGAACATGGCAGTGCAAGTG
Sequence-based reagentMus musculus Hoxd12 reverse primerThis paperPCR primersAGCTCTAGCTAGGCTCCTGTTTCATGC
Sequence-based reagentMus musculus Ptch1 forward primerThis paperPCR primersGGGAAGGCAGTTCATTGTTACTGTAACTG
Sequence-based reagentMus musculus Ptch1 reverse primerThis paperPCR primersTGTAATACGACTCACTATAGGTCAGAAGCTGCCACACACAGGCATGAAGC
Sequence-based reagentChiloscyllium punctatum Hand2 forward primerThis paperPCR primersACCAGCTACATTGCCTACCTCATGGAC
Sequence-based reagentChiloscyllium punctatum Hand2 reverse primerThis paperPCR primersCACTTGTTGAACGGAAGTGCACAAGTC
Sequence-based reagentChiloscyllium punctatum Vcan forward primerThis paperPCR primersAGCTTGGGAAGATGCAGAGAAGGAATG
Sequence-based reagentChiloscyllium punctatum Vcan reverse primerThis paperPCR primersAGAGCAGCTTCACAATGCAGTCTCTGG
Sequence-based reagentChiloscyllium punctatum Aldh1a2 forward primerThis paperPCR primersTTGAACTTGTACTAAGTGGTATCGCTG
Sequence-based reagentChiloscyllium punctatum Aldh1a2 reverse primerThis paperPCR primersAGGATGTGAACATTAGGCTGACCTCAC
Sequence-based reagentChiloscyllium punctatum Hoxd12 forward primerThis paperPCR primersGCCAGTATGCAACAGATCCTCTGATGG
Sequence-based reagentChiloscyllium punctatum Hoxd12 reverse primerThis paperPCR primersCTAATGACCTGTTGTACTTACATTCTC
Sequence-based reagentChiloscyllium punctatum Ptch1 forward primerThis paperPCR primersTTCAGCCAGATTGCAGATTACATCAACC
Sequence-based reagentChiloscyllium punctatum Ptch1 reverse primerThis paperPCR primersTTCTCTGTGTTTCACATTCAACGTCCTG
Commercial assay or kitNextera DNA Sample Preparation KitIlluminaCat# FC-121–1031
Commercial assay or kitTruSeq Stranded mRNA LT Sample Prep KitIlluminaCat# RS-122–2101
Software, algorithmTrinity
Software, algorithmBowtie2
Software, algorithmBWA
Software, algorithmMACS2
Software, algorithmHOMER
Software, algorithmRSEM
Software, algorithmscikit-learn

Additional files

Supplementary file 1

Summary of short-read sequencing data.
Supplementary file 2

Orthology asignment for the transcriptome of the brown-banded bamboo shark.

Columns 1–4: transcriptome assembly ID, NCBI gene ID, gene symbol, blast score.
Supplementary file 3

Orthology asignment for the gene model of the brown-banded bamboo shark.

Columns 1–4: gene model ID, NCBI gene ID, gene symbol, blast score.
Supplementary file 4

Quality control of orthology assignment.

Source data to create Figure 1C.
Supplementary file 5

The mean and SEM of TPM values of mouse limb RNA-seq data.

Source data for Figure 1D and other plots related to gene expression amount.
Supplementary file 6

The mean and SEM of TPM values of bamboo shark fin RNA-seq data.

Source data for Figure 1D and other plots related to gene expression amount.
Supplementary file 7

Clustered gene expression table with phenotype annotation.

The details of Figure 2A.
Supplementary file 8

The list of genes downregulated over time in mouse limb buds being upregulated in bamboo shark fin buds over time (related to Figure 2—figure supplement 1B).
Supplementary file 9

PCA loadings of Figure 4D and E.
Supplementary file 10

List of public data used in Figures 4, 5 and 6.
Supplementary file 11

GO analysis of ATAC-seq peaks.

c1 to c8 correspond to the clusters in Figure 5A.
Transparent reporting form

Download links

A two-part list of links to download the article, or parts of the article, in various formats.

Downloads (link to download the article as PDF)

Open citations (links to open the citations from this article in various online reference manager services)

Cite this article (links to download the citations from this article in formats compatible with various reference manager tools)

  1. Koh Onimaru
  2. Kaori Tatsumi
  3. Chiharu Tanegashima
  4. Mitsutaka Kadota
  5. Osamu Nishimura
  6. Shigehiro Kuraku
Developmental hourglass and heterochronic shifts in fin and limb development
eLife 10:e62865.