Skip to main content

Accelerated evolution and diversifying selection drove the adaptation of cetacean bone microstructure



The transition from land to sea by the ancestor of cetaceans approximately 50 million years ago was an incredible evolutionary event that led to a series of morphological, physiological, and behavioral adaptations. During this transition, bone microstructure evolved from the typical terrestrial form to the specialized structure found in modern cetaceans. While the bone microstructure of mammals has been documented before, investigations of its genetic basis lag behind. The increasing number of cetaceans with whole-genome sequences available may shed light on the mechanism underlying bone microstructure evolution as a result of land to water transitions.


Cetacean bone microstructure is consistent with their diverse ecological behaviors. Molecular evolution was assessed by correlating bone microstructure and gene substitution rates in terrestrial and aquatic species, and by detecting genes under positive selection along ancestral branches of cetaceans. We found that: 1) Genes involved in osteoclast function are under accelerated evolution in cetaceans, suggestive of important roles in bone remodeling during the adaptation to an aquatic environment; 2) Genes in the Wnt pathway critical for bone development and homeostasis show evidence of divergent evolution in cetaceans; 3) Several genes encoding bone collagens are under selective pressure in cetaceans.


Our results suggest that evolutionary pressures have shaped the bone microstructure of cetaceans, to facilitate life in diverse aquatic environments.


A bony skeleton is vital to many adaptive phenotypes in vertebrates and represented a major leap in evolution [1]. The microstructure and organization in bone reflects the biomechanical constraints that organisms undergo and generally show a strong ecological signal [2,3,4,5]. Bone microstructure has been used to infer the habitat and locomotor mode of extinct taxa, as well as to assess the ecological, biomechanical, and phylogenetic significance of bone microstructure of aquatic amniote groups [2, 4, 6,7,8]. Cetaceans (whales, dolphins, and porpoises) are the most speciose order of marine mammals (~ 89 extant species in 14 families) and inhabit diverse habitats, including ocean basins and large riverine ecosystems [9]. Extant cetaceans exhibit osteological adaptations to an aquatic lifestyle, accompanied with complex buoyancy control systems [10,11,12,13]. Taken together, cetaceans represent a model group for the study of bone microstructure adaptations by ecological transition.

Extensive anatomical records on bone microstructure exist. The link between bone microanatomy (e.g. limbs, vertebrae, and ribs) and habitat has been studied by various investigators [8, 10, 11, 13,14,15,16,17,18,19]. In general, flying taxa (e.g. bats) exhibit a ‘simple’ bone microanatomy, with thin cortices and few trabeculae in the medullary region [20]. Terrestrial mammals usually display an intermediate cortical thickness compared to aquatic mammals [19]. The specialized bone microstructure in most extant whale exhibits a thin layer of compact cortex, lacks a medulla, and has been described as ‘osteoporotic-like’ but without bone mass decrease or any pathological connotation [2]. Other types of bone microstructure, such as non-pathological densification (osteosclerosis) and swelling (pachyostosis) with increased bone density, are observed in relatively inactive shallow water dwellers such as Sirenia and Archaeocetes [8, 10, 21]. These specializations have been attributed to hydrodynamic or hydrostatic control of buoyancy, as well as diving and swimming [10, 22, 23].

The bony skeleton of vertebrates is a dynamic and metabolically living organ constituted primarily of calcium-phosphate minerals and type I collagen. The growth of bones is sculpted by modeling and continuously renewed by remodeling [1]. The bone remodeling process, which occurs throughout a lifetime to maintain mineral homeostasis, involves timed expression of osteoclasts and osteoblasts to balance the bone matrix [24]. Multiple genes and/or pathways must have been involved in the genetic adaption of bone microstructure in aquatic mammals. However, unlike the relatively well documented morphological changes of bone microstructure in cetaceans, the underlying genetic basis has not been well addressed in the literature. Multiple cetacean genomes are now available, allowing this question to be answered. To examine the evolution of bone microstructure in cetaceans, we contrasted the bone microstructure of cetaceans to other mammals and identified associated gene-phenotype correlations, and assessed the selective pressure of bone-remodeling genes on the ancestral lineage of cetaceans.


Sample collection and bone microstructure measurements

Anatomical data on bone microstructures (ribs 82 species, humeri of 14 species, and vertebrae of 50 species) (Additional file 1: Table S1–3) were obtained from two ways: (1) We collected and generated rib and humerus data from eight cetaceans and four adult terrestrial mammals in our lab – the false killer whale (Pseudorca crassidens), long-beaked common dolphin (Delphinus capensis), minke whale (Balaenoptera acutorostrata scammony), pantropical spotted dolphin (Stenella attenuate), Chinese white dolphin (Sousa chinensis), baiji (Lipotes vexillifer), finless porpoise (Neophocaena asiaeorientalis), common bottlenose dolphin (Tursiops truncatus), cow (Bos taurus), pig (Sus scrofa), sheep (Ovis aries), and dog (Canis lupus familiaris). Only dead stranding cetaceans for unknown reason in the wild was used in this study, and other mammal bone tissues were purchased in the market. No ethical approval was required. Bone sections were prepared as outlined by Canoville et al. and by Hayashi and colleagues [12, 19] (Fig. 1a). (2) To improve the sample size, sample images and information of more species was retrieved from previous various studies [2, 6, 12, 13, 19, 25].

Fig. 1

Characterization of bone microstructure. a Technical processing of the three kinds of bone, including long bone (rib and humerus) and vertebrae. Vertebrae data were obtained from [17], rib and humerus were sampled at the cross section at mid-length (solid line). Images of three bones were obtained from our laboratory. b Compactness profile indices for rib and humerus (Cg, Cc, Cp, S and P). Measured using Bone Profiler [25]. Each cross-section picture was converted to a binary image using Adobe Photoshop CS6

We employed binary images of thin sections and assessed bone density using Bone Profiler [25], to obtain the variables S (relative width of the transition zone between the medullary and the cortical regions) and P (proportional to the size of the medullary cavity) for each section. The compactness of the center/periphery/whole of bone sections was calculated using Image-Pro Plus, to obtain the variables Cg (globe compactness, ranging from 0 ~ 1), Cc (compactness in the center of the section. Cc values > 0 usually indicates the presence of trabeculae in the center of a section), and Cp (compactness in the periphery of a section. Cp values < 1 indicates that the cortical region displays porosity). All of these variables represent the proportion of mineralized bone tissue occupying the total sectional area (Fig. 1b). The parameter MD is the maximal diameter of bone section that can be used as a proxy for body size in statistical analysis of rib and humerus microanatomical data [6, 19].

Measurements of vertebrae from 50 species were retrieved from Houssaye and colleagues [17]. Measurements where more than half of the values were missing were excluded from our analysis. We included nine variables in the subsequent analyses: Cls (global compactness of the centrum in longitudinal section), CtsC (centrum compactness in transverse section), TNCL (total number of cavities in longitudinal section), NTCL (number of trabeculae in the centrum longitudinal section), AMCT (absolute mean cortical thickness in transverse section), RMCT (relative mean cortical thickness in transverse section), AMTT (absolute mean trabecular thickness), RMTT (relative mean trabecular thickness) and CL (centrum length used as proxy of body size in statistical analysis of vertebrae microanatomical data).

Statistical analysis of bone microanatomical data

Principal component analysis (PCA) was employed in order to reduce the dimensions of vertebrae indices. Vpc1 and Vpc2 (the first two principal components) explained 79% of the total variance (53 and 26%, respectively) of vertebrae compactness. Thus, Vpc1 and Vpc2 for vertebrae, and S, P, Cc, Cp, and Cg for rib variables were next used in phylogenetic generalized least squares (PGLS) multiple regression against body size variables (CL and MD for vertebrae and rib variables, respectively). Phylogenetic ANOVA analysis were employed to assess differences between habitat of bone variables for humerus, rib, and vertebrae using the ‘phytools’ package in R [26].

For the rib data set (n = 82), we further divided aquatic habitat species into shallow water/coastal water swimmers, and deep divers based on ecological behavior characteristics and diving depth data [27]. To assess bone histological parameters of 24 marine mammals we compared the S/P/Cc/Cp/Cg of ribs in six habitat categories in box plot and performed PGLS regression analysis of these variables against diving depth (in meters).

Orthologous preparation and phylogenetic comparative analyses

Of the species used in the quantitative analysis of bone microanatomy, 27 (or a related species in the same family or similar habitat or locomotion) had corresponding whole-genome sequences (Fig. 2). We generated three bone microstructure data sets: ribs for 27 species, humeri for 14 species, and vertebrae for 13 species (Additional file 1: Table S4). We next obtained genome sequences from NCBI and obtained 1:1 orthologous genes among them using OrthoMCL v2.0.9 [28]. Of 3621 single-copy orthologs, 348 were classified as bone function-associated genes based on GO terms (‘skeletal system development’, ‘ossification’, ‘bone remodeling’, ‘osteoblast proliferation’, ‘osteoclast differentiation’, and ‘osteoclast proliferation’), KEGG pathway names (‘osteoclast differentiation’), and a literature survey (key words ‘bone development’ and/or ‘osteoclast’) [29, 30]. Multiple alignments of orthologous sequences were generated using PRANK v150803 [31], followed by Gblocks v0.91b [32] and manual curation of alignments.

Fig. 2

Phylogenetic tree of 27 species used in regression analyses. Representatives from marine, terrestrial, and flying mammals are in blue, black, and green, respectively. Binary cross section images of ribs with variation in five compactness indices were shown for each species. Bone section images were obtained from [19] and our lab: Tursiops truncatus (PL15–0145-28), Orcinus orca (AMNH 34261), Neophocaena asiaeorientalis (PL15–0145-21), Lipotes vexillifer (PL15–0145-18), Balaenoptera acutorostrata (PL15–0145-08), Bison bison (MHNL 50002450), Bos taurus (PL15–0145-33), Ovis aries (PL15–0145-39), Camelus bactrianus (MHNL 50002066), Camelus dromedaries (MHNL 50002063), Ceratotherium sinum (AMNH 51855), Equus caballus (MHNL 50002029); Felis catus (represent by Felis silvestris p.c. VB), Canis familiaris (PL15–0145-43), Ursus maritimus (p.c. VB), Odobenus rosmarus (MHNL 50001014), Leptonychotes weddellii (represent by Cystophora cristata AMNH 184659), Pteropus Alecto (represent by Eidolon helvum ZFMK no sp. number), Pteropus vampyrus (represent by Pteropus giganteus ZFMK 80.851), Eptesicus fuscus (represent by Rousettus aegyptiacus ZFMK 2001.004), Myotis davidii (represent by Pipistrellus pipistrellus ZFMK a), Erinaceus europaeus (p.c. VB), Gorilla gorilla gorilla (MHNL 50001762), Callithrix jacchus (ZFMK MAM_1983.0366), Rattus norvegicus (p.c. VB), Trichechus manatus (represent by Dugong dugon MHNL 50002521), Loxodonta africana (represent by Elephas maximus MHNL 50002671)

The ‘root-to-tip’ dN/dS, defined as the average value of accumulated dN/dS extending from the last common ancestor of all mammals examined to the respective terminal branch, was estimated with PAML v4.4 [33] using a free-ratio model and parsed using custom Perl scripts. This measurement has been recognized as an index of selection which takes the entire evolutionary history of a lineage from a common ancestor into account and negates the issue of temporal effects on dN/dS [34]. Regression between ‘root-to-tip’ dN/dS and bone histological parameters was assessed using PGLS regression models under a phylogenetic framework across mammals (using the R package ‘caper’ v0.5.2) [35]. Topology and divergence date of a 27-species phylogenetic tree was obtained from the online resource TimeTree [36]. Briefly, PGLS employs a phylogenetic tree as input to assess the impact of phylogenetic non-independence between species. For every regression analysis, a quantitative measure of phylogenetic signal (Pagel’s lambda; λ) is calculated through maximum likelihood estimations. A λ value of 1, of or near 1, indicates that a variable is fully explained by evolutionary history and thus shows a strong phylogenetic signal [37, 38]. To obtain more stringent correlation P values, we further employed a two-step verification procedure [39]. On the basis of ‘P value.all’ from the regression analyses for all samples, the following two P values were calculated: 1) ‘P value.robust’ from PGLS repeated after excluding the sample with largest residual error; 2) ‘P value.max’ from the regression on the remaining species, to calculate the maximal P value after dropping one species.

Estimating selective pressure

To identify genes that had been under selection, we set the ancestral branch of cetaceans as the foreground branch. Two models, the branch-site model [40] and clade model C [41], were implemented using codeml in PAML v4.4 [33]. The P value of each gene was computed using likelihood ratio tests (LRTs). The clade model C can detect evidence of divergent selective pressures acting across the cetacean clade as the foreground compared with the remaining species in the tree as the background. We set each model with three initial ω values (0.5, 1, and 1.5), to obtain the robust average ω, and compared this result with model M2a_ref (nearly neutral) via LRTs. Only genes with a unchangeable likelihood value for three initial ω values were considered interesting. To identify associations between genes, we implemented STRING (v10.5) functional analysis [42]. STRING integrates predicted and experimentally confirmed relationships between proteins that are likely to contribute to a common biological purpose.


Ecological signals of bone microstructure across mammals

For three anatomical datasets (Additional file 1: Table S1–3), significant positive correlations were found between body size and the rib variables Cc (P = 3.9 × 10− 9, r2 = 0.35), Cp (P = 1.03 × 10− 3, r2 = 0.16), S (P = 2.15 × 10− 10, r2 = 0.39), and between body size and the vertebrae variable Vpc1 (P = 2.2 × 10− 16, r2 = 0.76) (Additional file 2: Figure S1). We classified habitats into five categories (aquatic, amphibious, arboreal, terrestrial, and flying. For the humerus, only aquatic and terrestrial were considered). Phylogenetic ANOVA analysis showed that the aquatic group was significantly different from other groups in terms of S, Cc, and Cp (P < 0.01 in rib; P < 0.05 in humerus) (Additional file 2: Table S5).

Of the 24 marine mammals examined in rib data set (Additional file 1: Table S1), deep divers (false killer whale, Pseudorca crassidens; narwhal, Monodon monoceros; Blainville’s beaked whale, Mesoplodon densirostris; long-finned pilot whale, Globicephala melas; and hooded seal, Cystophora cristata) had lower Cp, Cc, S and Cg, and higher P (Additional file 2: Figure S2 and S3). Additionally, diving depth (max/average) significantly correlated with Cp (max depth; P < 0.05, r2 = 0.13; average depth; P < 0.05, r2 = 0.13) and Cg (max depth; P < 0.05, r2 = 0.13; average depth; P < 0.05, r2 = 0.19) (Additional file 2: Figure S4), suggesting that deep-diving species have a lower bone density. For example, the hooded seal (C. cristata) and Blainville’s beaked whale (M. densirostris) exhibited low global compactness and numerous thin bone trabeculae.

Detection of selective pressure and gene–phenotype association analysis

The branch-site model was used to identify positively selected genes in a 27-species data set with genome and corresponding bone microstructure measurements (see Additional file 1: Table S4 and Fig. 2). Nine genes (COL1A2, COL3A1, FSHR, IFNAR1, MEPE, MITF, NFATC3, TEC, and TNFRSF1A) showed evidence of strong positive selection (likelihood ratio tests, LRTs P < 0.05) in the common ancestor of cetaceans. A total of 14 genes (SPARC, COL2A1, COL9A1, COL5A2, COL3A1, HES1, CTNNB1, FZD4, RUNX2, DVL3, TEC, FOSL1, STAT1, and LCP2) showed significant positive selection in cetaceans but not in the outgroup taxa (foreground ω2 > 1.00 and background ω1 < 1.00, LRT P < 0.05). Interestingly, assessment for divergent evolution (clade model C) identified genes related to ‘Wnt signaling’ and ‘regulation of osteoblasts’ (Table 1 and see Additional file 3: Table S6–7).

Table 1 Summary of genes with a ‘root-to-tip’ dN/dS significantly correlated with indices of three kinds of bone (rib, humerus and vertebrae) and under positive selection (PSG) or divergent selection (DSG) in cetaceans

A comparison of ‘root-to-tip’ dN/dS values (see Additional file 1: Table S4) and bone measurements from the ribs (27 species), humeri (14 species), and vertebrae (13 species) revealed 83 genes with substitution rates correlating with bone variables (see Additional file 4: Table S8–10). Since body size correlated positively with several variables (S, Cc, and Cp for rib; Vpc1 for vertebrae) (Additional file 2: Figure S1), variable residuals were computed to account for body size influence on regression and employed in a new round of regression analysis. For the variables Cc and Cp, 83% of genes were significant or near significant (P value.all for residual < 0.1), while one gene was detected near significant each for the variables S (S ~ NCF2: P value.all for residual = 0.08) and Vpc1 (vpc1~SYK: P value.all for residual = 0.083) (Additional file 4: Table S8). Of the 83 genes, 27 genes are associated with the KEGG pathway ‘Osteoclast differentiation’ (P < 0.1) (Fig. 3) and six were with ‘Phagosome’ (P < 0.01) (see Additional file 5: Table S11–12). Moreover, six genes (PIK3CB, NGF, MAPK13, TNFRSF11B, NCF2, and IFNGR2) associated with the osteoclast differentiation pathway showed a significant correlation with two out of three bone types (Fig. 4a). Four genes (MITF, TNFRSF1A, TEC, and IFNAR1) in the osteoclast differentiation pathway were identified by the branch site model, suggesting that they have undergone positive selection in cetaceans. Additionally, clade model C revealed that four genes (STAT1, TEC, LCP2, and FOSL1) have been subject to divergent selection pressures in cetaceans compared to other taxa (Fig. 3 and 4a; Table 1). With the exception of TNFRSF1A (TNF receptor superfamily member 1A), 26 genes enriched for the osteoclast differentiation pathway showed a negative correlation between ‘root-to-tip’ dN/dS and Cp and a positive correlation between ‘root-to-tip’ dN/dS and bone indices Cc and S, (Fig. 4b). Thus, genes associated with osteoclast function have evolved rapidly in aquatic mammals.

Fig. 3

Mapping of genes with ‘root-to-tip’ dN/dS significantly correlated with bone microstructures or under positive and divergent selections to the KEGG osteoclast differentiation pathway

Fig. 4

Overview of 83 genes with gene substitution rates correlating with bone microstructures. a. Protein–protein interaction network generated using STRING [42]. Nodes represent protein-coding genes correlating with bone variables (regression analysis, vs. ‘root-to-tip’ dN/dS). Lines between nodes indicate inferred or experimentally demonstrated biological associations. Humerus, rib, and vertebrae are indicated in yellow, beige, and brown respectively. Genes also under positive or divergent selection are indicated in maroon and purple, respectively. b. Heat map of genes correlating (green, positive; red, negative) with bone variables (regression analysis, vs. ‘root-to-tip’ dN/dS). Genes in the osteoclast differentiation pathway are highlighted in pink. For selected genes, a plot of bone variable (vertical axis; please see main text for details) and substitution rate (‘root-to-tip’ dN/dS) is shown; with salmon, blue, and green indicating species in an aquatic, flying, and terrestrial habitat, respectively

The association of ‘root-to-tip’ dN/dS with bone mass parameters and the selective pressure analysis both identified several genes which encode collagen proteins (Table 1 and Fig. 4a), essential components of bone and cartilage and required for endochondral ossification. For example, genes encoding type I and III collagens (COL1A2 and COL3A1) were under positive selection in the ancestral lineage of cetaceans. Although our two-step verification procedure (see [39]) revealed a relationship between S and Cc vs dN/dS of COL1A2 (Cc: P value.all = 0.002, P value.robust = 0.001, P value.max = 0.002; S: P value.all = 0.004, P value.robust = 0.007, P value.max = 0.049), this was not significant in regression analyses where variable residuals were employed to account for body size (Cc: P = 0.305; S: P = 0.238). Furthermore, the divergent selection scan showed that type II, V, and IX collagens (COL2A1, COL5A2 and COL9A1) and osteonectin (SPARC) were also under positive selection in extant cetaceans.


Bone microstructure in cetaceans with diverse ecological behaviors

Anatomical records from marine mammals (16 cetaceans, seven pinnipeds, and the sea otter) allowed us to address the great diversity of bone microstructure in marine mammals, especially for the fully aquatic cetaceans. Combining qualitative data and statistical analyses, we reveal that whole bone volume (Cg) does not scale substantially with habitat or body size. Mammals with a larger body size tend to have a more complex organization, with an increased transition zone (S) between the cortical bone and the medullary cavity, as well as numerous bone trabeculae occupying the medullary region (also called cancellous bone; Cc). For animals from different categories of habitats, the compactness of cortical bone (Cp) is relatively low in aquatic mammals, especially in deep-diving whales. Cortical bone is dense and strong enough to provide both support and protection for most mammals. However, for cetaceans, who harbor a more porous bone structure, it seems that the trabecular bone tends to invade into the cortical bone – as evidenced by higher S and Cc values. It is generally accepted that the trabecular bone is more responsive to, and malleable for, variations in magnitude and direction of load throughout life compared to cortical bone [59]. Therefore, this pattern of bone microstructure may be an adaptation for the aquatic lifestyle of cetaceans.

The inner bone structure within cetaceans is variable, especially in the case of the ribs. It has been assumed that ribs contribute significantly to the mass and inertia of the body in these “limbless” mammals [19]. Compared to most terrestrial mammals, which have a fully compact cortical bone, cetaceans display porosity (lower Cp). Some cetaceans in the Delphinidae and Monodontidae families exhibit unexpectedly thick cortices, with a certain extent of porosity and thicker bone trabeculae in the medullary region. This could be because an increased bone density supports feats such as the acrobatics and fast swimming speed of dolphins and feeding in deep divers [60, 61]. Another interesting observation is that species which reside in shallow freshwater, such as baiji (Lipotes vexillifer) and boto (Inia geoffrensis), have a spongy section with numerous thin bone trabeculae and a low global compactness (Cg). Although the function of this modification is unknown, it may support relatively lower buoyancy in a freshwater habitat.

Accelerated evolution of osteoclast differentiation-related genes in cetaceans

Previous studies suggested that down- and up-regulation of osteoclast activity was central to bone microstructure of tetrapods returning to an aquatic environment [10]. Our results reveal a positive correlation to S or Cc and negative correlation to Cp by genes in the osteoclast differentiation pathway. Since aquatic mammals had a relatively higher S and Cc, and a lower Cp, genes with higher ‘root-to-tip’ dN/dS values in cetaceans could reflect accelerated evolution of osteoclastogenesis.

Most of the positively selected genes (e.g., TEC, LCP2, TRAF2, MITF, CTSK, LCK, and GRB2) in the ancestral branch of cetaceans are genes with ‘root-to-tip’ dN/dS that negatively correlate with Cp. These genes encode intracellular signaling cascade proteins involved in osteoclast differentiation (see Fig. 4b and Table 1). Mutations in, or the lack of any of these, genes would cause severe osteopetrosis [54, 57, 62,63,64]. The genes (IFNAR1, TNFRSF1A, TNFRSF11B, TNFSF11, and TGFB2) in the osteoclast differentiation pathway are ligands or receptors that trigger the differentiation process of osteoclasts (Fig. 3 and Table 1). Two genes have a negative impact on osteoclast differentiation. The ‘root-to-tip’ dN/dS of TNFRSF11B (osteoprotegerin; also known as OPG) positively correlated with Cc in both rib (P value.robust < 3.89 × 10− 06, r2 = 0.91) and humerus (P value.robust < 2.10 × 10− 06, r2 = 0.6), while TNFSF11 (receptor activator of nuclear factor kappa B ligand; also known as RANKL) negatively correlated to Vpc2 in vertebrae (P value.robust < 1.61 × 10− 04, r2 = 0.75). In GWAS meta-analyses these genes were associated with volumetric bone mineral density, suggesting that the RANK-RANKL-OPG axis affects the skeleton – at least in part by influencing the density of cortical bone [65]. Another positively selected gene in the cetacean lineage of interest is TNFRSF1A (also known as TNFR1), the only gene with a ‘root-to-tip’ dN/dS negatively correlated with Cc (P value.robust = 3.93 × 10− 05, r2 = 0.49). TNFRSF1A is the receptor of TNF-α and modulates immune and inflammatory processes, as well as bone homeostasis [58]. Cross-talk between interferons and other cytokines in bone remodeling have recently received greater attention. For example, one of the IFN-β receptor components in the type I interferon system is IFNAR1 (Table 1). Ifnar1−/− mice have markedly reduced trabecular bone mass, a key feature of osteoporosis, suggesting that the gene plays a critical role in osteoclastogenesis [56]. We hypothesize that bone homeostasis in aquatic mammals is maintained by the RANK-RANKL-OPG axis acting in concert with cytokines.

In summary, we have identified a correlation between accelerated changes of osteoclast-associated genes and increased trabeculae in cetaceans. This correlation could reflect increased bone resorption. Further analyses of selective pressure suggest that several genes related to osteoclast differentiation underwent accelerated change due to positive selection, likely to allow bone microstructure specialization in cetaceans.

Divergent selection of bone formation genes between cetaceans and other mammals

Genes associated with ‘canonical Wnt signaling’ and ‘regulation of osteoblasts’ were shown to evolve through divergent selection in cetaceans compared to other taxa (Table 1). It is well known that Wnt signaling plays a pivotal role in skeletal homeostasis by regulating bone formation and bone resorption by osteoblasts and osteoclasts [43]. Genes associated with the Wnt intracellular signaling (CTNNB1, FZD4, and DVL3) were under divergent selection in cetaceans. Other genes under divergent selection encode transcription factors. This includes RUNX2 which is essential for the maturation of osteoblasts [50]. Knockout of RUNX2 in mice results in a complete lack of bone formation and arrested osteoblast differentiation [66]. Moreover, studies of mouse models with knockout or transgenic Wnt pathway components have demonstrated that this signaling pathway regulates most aspects of osteoblast physiology, including bone matrix formation/mineralization, osteoclastogenesis, and bone resorption [43, 44]. Thus, significant divergent selection of these genes between cetaceans and other mammals suggests that they played crucial roles in driving the bone development and formation in response to aquatic adaptations.

Positive selection of collagen genes in cetaceans

The matrix of bone mainly contains collagen fibers and mineral deposits. Type I collagen in bone and type II collagen in cartilage provide structural integrity and account for mechanical strength [67]. Osteogenesis imperfecta, a severe genetic disorder manifested by increased bone fragility and low bone mass, is included in the diseases associated with COL1A2 [46]. Alteration of collagen structural properties by missense mutations in COL1A1 and COL1A2 reduce the mechanical properties of bone [68]. Moreover, cetacean COL1A2 harbors unique substitutions that may modify the collagen triple helix [69]. Nevertheless, the fact that COL1A2 and several other genes mentioned before were not statistical significant after correcting for body size residuals. Possible reasons may include sample size or ecological/statistical implications for bone sample feature, e.g. only one specific region of the rib series were sampled or intraspecific microanatomical variability, which were mentioned in Canoville et al. (2016) [19]. Type III collagen (COL3A1) showed a negative correlation between evolutionary rate and Cc in cetaceans, as well as positive selection and divergent selection in the ancestral lineage of cetaceans. COL3A1 is highly expressed during embryonic skeletal development and expressed by osteoblasts in mature bone. Both in vivo and in vitro experiments have shown that COL3A1 plays an important role in the development of trabecular bone through its effects on osteoblast differentiation [70]. Loss of COL3A1 function is associated with Ehlers-Danlos syndrome (EDS), where skeletal manifestations include a distinctive facial appearance, scoliosis, and osteoporosis [71]. Cartilaginous collagen fibrils are represented by collagens type II, IX, and XI and endochondral ossification is gradually replaced by a bone matrix composed of type I and type II collagens [72]. Together with divergently selected osteonectin (SPARC), a gene which encodes a protein that bind to ECM proteins such as type I, III, IV, and V collagens to promote bone mineralization [47, 73], we speculate that bone collagen are likely to have contributed significantly to the specialization of cetaceans bone microstructure.


Along with the transition of cetaceans from land to sea, changes in habitats and ecological behaviors in cetaceans resulted in diverse patterns of bone microstructures. In the present study, we found that genes involved in osteoclast differentiation experienced accelerated evolution in cetaceans. In addition, Wnt signaling and osteoblastogenesis associated with bone development were found to be under divergent evolution in cetaceans. We also detected positive selection of collagen proteins in cetaceans. These findings provide new insights into the nexus between genes and the secondary adaptation of terrestrial mammals to aquatic life.

Availability of data and materials

The data generated and analyzed during this study are included in this article and its additional files 1-5.



American Museum of Natural History, New York, NY, USA


Gene Ontology


uncatalogued research collection of John R. Hutchinson at The Royal Veterinary College, Hatfield, UK


Musée des Confluences, Centre de conservation et d’étude des collections, 13A rue. Bancel, 69,007 Lyon, France


Museum National d’Histoire Naturelle, Paris, France


ordinary least squares regression models

p.c. VB:

Personal collection of V. de Buffrénil


phylogenetic generalized least squares. Institutional abbreviations


Zoological Research Museum Alexander Koenig, Bonn, Germany


  1. 1.

    Martin RB, Burr DB, Sharkey NA, Fyhrie DP. Skeletal tissue mechanics: springer; 2015.

    Google Scholar 

  2. 2.

    Houssaye A, Martin Sander P, Klein N. Adaptive patterns in aquatic amniote bone microanatomy—more complex than previously thought. Integrative comparative biology. 2016;56(6):1349–69.

    PubMed  Article  Google Scholar 

  3. 3.

    Ruimerman R, Van Rietbergen B, Hilbers P, Huiskes R. The effects of trabecular-bone loading variables on the surface signaling potential for bone remodeling and adaptation. Ann Biomed Eng. 2005;33(1):71–8.

    CAS  PubMed  Article  Google Scholar 

  4. 4.

    Cooper LN, Clementz MT, Usip S, Bajpai S, Hussain ST, Hieronymus TL. Aquatic habits of cetacean ancestors: integrating bone microanatomy and stable isotopes. Integrative comparative biology. 2016;56(6):1370–84.

    PubMed  Article  Google Scholar 

  5. 5.

    Chappard D, Baslé M-F, Legrand E, Audran M. Trabecular bone microarchitecture: a review. Morphologie. 2008;92(299):162–70.

    CAS  PubMed  Article  Google Scholar 

  6. 6.

    Houssaye A, Waskow K, Hayashi S, Cornette R, Lee AH, Hutchinson JR. Biomechanical evolution of solid bones in large animals: a microanatomical investigation. Biol J Linn Soc. 2016;117(2):350–71.

    Article  Google Scholar 

  7. 7.

    Amson E, de Muizon C, Laurin M, Argot C, de Buffrénil V. Gradual adaptation of bone structure to aquatic lifestyle in extinct sloths from Peru. Proc R Soc Lond B Biol Sci. 2014;281(1782):20140192.

    Article  Google Scholar 

  8. 8.

    Houssaye A, Tafforeau P, De Muizon C, Gingerich PD. Transition of Eocene whales from land to sea: evidence from bone microstructure. PLoS One. 2015;10(2):e0118409.

    PubMed  PubMed Central  Article  CAS  Google Scholar 

  9. 9.

    Jefferson TA, Webber MA, Pitman RL. Marine mammals of the world: a comprehensive guide to their identification: Elsevier; 2011.

    Google Scholar 

  10. 10.

    Gray NM, Kainec K, Madar S, Tomko L, Wolfe S. Sink or swim? Bone density as a mechanism for buoyancy control in early cetaceans. The Anatomical Record: Advances in Integrative Anatomy Evolutionary Biology: Advances in Integrative Anatomy Evolutionary Biology. 2007;290(6):638–53.

    Article  Google Scholar 

  11. 11.

    Dumont M, Laurin M, Jacques F, Pelle E, Dabin W, de Buffrénil V. Inner architecture of vertebral centra in terrestrial and aquatic mammals: a two-dimensional comparative study. J Morphol. 2013;274(5):570–84.

    PubMed  Article  Google Scholar 

  12. 12.

    Hayashi S, Houssaye A, Nakajima Y, Chiba K, Ando T, Sawamura H, Inuzuka N, Kaneko N, Osaki T. Bone inner structure suggests increasing aquatic adaptations in Desmostylia (Mammalia, Afrotheria). PLoS One. 2013;8(4):e59146.

    CAS  PubMed  PubMed Central  Article  Google Scholar 

  13. 13.

    Laurin M, Canoville A, Germain D. Bone microanatomy and lifestyle: a descriptive approach. Comptes Rendus Palevol. 2011;10(5–6):381–402.

    Article  Google Scholar 

  14. 14.

    Canoville A, Laurin M. Microanatomical diversity of the humerus and lifestyle in lissamphibians. Acta Zool. 2009;90(2):110–22.

    Article  Google Scholar 

  15. 15.

    Canoville A, Laurin M. Evolution of humeral microanatomy and lifestyle in amniotes, and some comments on palaeobiological inferences. Biol J Linn Soc. 2010;100(2):384–406.

    Article  Google Scholar 

  16. 16.

    Quemeneur S, De Buffrenil V, Laurin M. Microanatomy of the amniote femur and inference of lifestyle in limbed vertebrates. Biol J Linn Soc. 2013;109(3):644–55.

    Article  Google Scholar 

  17. 17.

    Houssaye A, Tafforeau P, Herrel A. Amniote vertebral microanatomy–what are the major trends? Biol J Linn Soc. 2014;112(4):735–46.

    Article  Google Scholar 

  18. 18.

    Waskow K, Sander PM. Growth record and histological variation in the dorsal ribs of Camarasaurus sp.(Sauropoda). J Vertebr Paleontol. 2014;34(4):852–69.

    Article  Google Scholar 

  19. 19.

    Canoville A, de Buffrénil V, Laurin M. Microanatomical diversity of amniote ribs: an exploratory quantitative study. Biol J Linn Soc. 2016;118(4):706–33.

    Article  Google Scholar 

  20. 20.

    Swartz SM, Bennett MB, Carrier DR. Wing bone stresses in free flying bats and the evolution of skeletal design for flight. Nature. 1992;359(6397):726.

    CAS  PubMed  Article  Google Scholar 

  21. 21.

    de Buffrénil V, Canoville A, D’Anastasio R, Domning DP. Evolution of sirenian pachyosteosclerosis, a model-case for the study of bone structure in aquatic tetrapods. J Mamm Evol. 2010;17(2):101–20.

    Article  Google Scholar 

  22. 22.

    Houssaye A. “Pachyostosis” in aquatic amniotes: a review. Integrative Zoology. 2009;4(4):325–40.

    PubMed  Article  Google Scholar 

  23. 23.

    George JC, Stimmelmayr R, Suydam R, Usip S, Givens G, Sformo T, Thewissen J. Severe bone loss as part of the life history strategy of bowhead whales. PLoS One. 2016;11(6):e0156753.

    PubMed  PubMed Central  Article  CAS  Google Scholar 

  24. 24.

    Doherty AH, Ghalambor CK, Donahue SW. Evolutionary physiology of bone: bone metabolism in changing environments. Physiology. 2015;30(1):17–29.

    CAS  PubMed  Article  Google Scholar 

  25. 25.

    Girondot M, Laurin M. Bone profiler: a tool to quantify, model, and statistically compare bone-section compactness profiles. J Vertebr Paleontol. 2003;23(2):458–61.

    Article  Google Scholar 

  26. 26.

    Revell LJ. Phytools: an R package for phylogenetic comparative biology (and other things). Methods in Ecology Evolution. 2012;3(2):217–23.

    Article  Google Scholar 

  27. 27.

    Watwood SL, Buonantony DM: Dive distribution and group size parameters for marine species occurring in navy training and testing areas in the north Atlantic and north Pacific oceans. In.: NAVAL UNDERSEA WARFARE CENTER DIV NEWPORT RI; 2012.

  28. 28.

    Li L, Stoeckert CJ, Roos DS. OrthoMCL: identification of ortholog groups for eukaryotic genomes. Genome Res. 2003;13(9):2178–89.

    CAS  PubMed  PubMed Central  Article  Google Scholar 

  29. 29.

    Kanehisa M, Sato Y, Kawashima M, Furumichi M, Tanabe M. KEGG as a reference resource for gene and protein annotation. Nucleic Acids Res. 2015;44(D1):D457–62.

    PubMed  PubMed Central  Article  CAS  Google Scholar 

  30. 30.

    Consortium GO: The Gene Ontology (GO) database and informatics resource. Nucleic acids research 2004, 32(suppl_1):D258-D261.

  31. 31.

    Löytynoja A, Goldman N. Phylogeny-aware gap placement prevents errors in sequence alignment and evolutionary analysis. Science. 2008;320(5883):1632–5.

    PubMed  Article  CAS  Google Scholar 

  32. 32.

    Castresana J. Selection of conserved blocks from multiple alignments for their use in phylogenetic analysis. Mol Biol Evol. 2000;17(4):540–52.

    CAS  Article  Google Scholar 

  33. 33.

    Yang Z. PAML 4: phylogenetic analysis by maximum likelihood. Mol Biol Evol. 2007;24(8):1586–91.

    CAS  PubMed  Article  Google Scholar 

  34. 34.

    Montgomery SH, Capellini I, Venditti C, Barton RA, Mundy NI. Adaptive evolution of four microcephaly genes and the evolution of brain size in anthropoid primates. Mol Biol Evol. 2010;28(1):625–38.

    PubMed  Article  CAS  Google Scholar 

  35. 35.

    Orme D. The caper package: comparative analysis of phylogenetics and evolution in R. R package version. 2013;5(2):1–36.

    Google Scholar 

  36. 36.

    Hedges SB, Dudley J, Kumar S. TimeTree: a public knowledge-base of divergence times among organisms. Bioinformatics. 2006;22(23):2971–2.

    CAS  PubMed  PubMed Central  Article  Google Scholar 

  37. 37.

    Symonds MR, Blomberg SP: A primer on phylogenetic generalised least squares. In: Modern phylogenetic comparative methods and their application in evolutionary biology. Springer; 2014: 105–130.

  38. 38.

    Pagel M. Inferring the historical patterns of biological evolution. Nature. 1999;401(6756):877.

    CAS  Article  Google Scholar 

  39. 39.

    Ma S, Yim SH, Lee S-G, Kim EB, Lee S-R, Chang K-T, Buffenstein R, Lewis KN, Park TJ, Miller RA. Organization of the mammalian metabolome according to organ function, lineage specialization, and longevity. Cell Metab. 2015;22(2):332–43.

    CAS  PubMed  PubMed Central  Article  Google Scholar 

  40. 40.

    Zhang J, Nielsen R, Yang Z. Evaluation of an improved branch-site likelihood method for detecting positive selection at the molecular level. Mol Biol Evol. 2005;22(12):2472–9.

    CAS  PubMed  Article  Google Scholar 

  41. 41.

    Weadick CJ, Chang BS. An improved likelihood ratio test for detecting site-specific functional divergence among clades of protein-coding genes. Mol Biol Evol. 2011;29(5):1297–300.

    PubMed  Article  CAS  Google Scholar 

  42. 42.

    Szklarczyk D, Morris JH, Cook H, Kuhn M, Wyder S, Simonovic M, Santos A, Doncheva NT, Roth A, Bork P. The STRING database in 2017: quality-controlled protein–protein association networks, made broadly accessible. Nucleic acids research. 2016:gkw937.

  43. 43.

    Glass DA, Bialek P, Ahn JD, Starbuck M, Patel MS, Clevers H, Taketo MM, Long F, McMahon AP, Lang RA. Canonical Wnt signaling in differentiated osteoblasts controls osteoclast differentiation. Dev Cell. 2005;8(5):751–64.

    CAS  PubMed  Article  Google Scholar 

  44. 44.

    Yavropoulou MP, Yovos JG. The role of the Wnt signaling pathway in osteoblast commitment and differentiation. HORMONES-ATHENS. 2007;6(4):279.

    PubMed  Article  Google Scholar 

  45. 45.

    Sun L, Peng Y, Sharrow AC, Iqbal J, Zhang Z, Papachristou DJ, Zaidi S, Zhu L-L, Yaroslavskiy BB, Zhou H. FSH directly regulates bone mass. Cell. 2006;125(2):247–60.

    CAS  PubMed  PubMed Central  Article  Google Scholar 

  46. 46.

    Fratzl P: Collagen: structure and mechanics, an introduction. In: Collagen. Springer; 2008: 1–13.

  47. 47.

    Kelm R, Mann K. The collagen binding specificity of bone and platelet osteonectin is related to differences in glycosylation. J Biol Chem. 1991;266(15):9632–9.

    CAS  PubMed  Google Scholar 

  48. 48.

    Ribeiro N, Sousa SR, Brekken RA, Monteiro FJ. Role of SPARC in bone remodeling and cancer-related bone metastasis. J Cell Biochem. 2014;115(1):17–26.

    CAS  PubMed  Article  Google Scholar 

  49. 49.

    Zanotti S, Smerdel-Ramoya A, Canalis E. HES1 (hairy and enhancer of split 1) is a determinant of bone mass. J Biol Chem. 2011;286(4):2648–57.

    CAS  PubMed  Article  Google Scholar 

  50. 50.

    Komori T. Regulation of bone development and maintenance by Runx2. Front Biosci. 2008;13:898–903.

    CAS  PubMed  Article  Google Scholar 

  51. 51.

    Komori T. Regulation of bone development and extracellular matrix protein genes by RUNX2. Cell Tissue Res. 2010;339(1):189.

    CAS  Article  Google Scholar 

  52. 52.

    Kulkarni RN, Bakker AD, Everts V, Klein-Nulend J. Inhibition of osteoclastogenesis by mechanically loaded osteocytes: involvement of MEPE. Calcif Tissue Int. 2010;87(5):461–8.

    CAS  PubMed  PubMed Central  Article  Google Scholar 

  53. 53.

    Kim S, Koga T, Isobe M, Kern BE, Yokochi T, Chin YE, Karsenty G, Taniguchi T, Takayanagi H. Stat1 functions as a cytoplasmic attenuator of Runx2 in the transcriptional program of osteoblast differentiation. Genes Dev. 2003;17(16):1979–91.

    CAS  PubMed  PubMed Central  Article  Google Scholar 

  54. 54.

    Shinohara M, Koga T, Okamoto K, Sakaguchi S, Arai K, Yasuda H, Takai T, Kodama T, Morio T, Geha RS. Tyrosine kinases Btk and Tec regulate osteoclast differentiation by linking RANK and ITAM signals. Cell. 2008;132(5):794–806.

    CAS  PubMed  Article  Google Scholar 

  55. 55.

    Matsuo K, Owens JM, Tonko M, Elliott C, Chambers TJ, Wagner EF. Fosl1 is a transcriptional target of c-Fos during osteoclast differentiation. Nat Genet. 2000;24(2):184.

    CAS  PubMed  Article  Google Scholar 

  56. 56.

    Takayanagi H, Kim S, Matsuo K, Suzuki H, Suzuki T, Sato K, Yokochi T, Oda H, Nakamura K, Ida N. RANKL maintains bone homeostasis through c-Fos-dependent induction of interferon-β. Nature. 2002;416(6882):744.

    CAS  PubMed  Article  Google Scholar 

  57. 57.

    Sharma SM, Bronisz A, Hu R, Patel K, Mansky KC, Sif S, Ostrowski MC. MITF and PU. 1 recruit p38 MAPK and NFATc1 to target genes during osteoclast differentiation. Journal of Biological Chemistry. 2007.

  58. 58.

    Osta B, Benedetti G, Miossec P. Classical and paradoxical effects of TNF-α on bone homeostasis. Front Immunol. 2014;5:48.

    PubMed  PubMed Central  Google Scholar 

  59. 59.

    Kivell TL. A review of trabecular bone functional adaptation: what have we learned from trabecular analyses in extant hominoids and what can we apply to fossils? J Anat. 2016;228(4):569–94.

    PubMed  PubMed Central  Article  Google Scholar 

  60. 60.

    Currey JD, Alexander RM. The thickness of the walls of tubular bones. J Zool. 1985;206(4):453–68.

    Article  Google Scholar 

  61. 61.

    Webb PW, De Buffrénil V. Locomotion in the biology of large aquatic vertebrates. Trans Am Fish Soc. 1990;119(4):629–41.

    Article  Google Scholar 

  62. 62.

    Kanazawa K, Kudo A. TRAF2 is essential for TNF-α-induced osteoclastogenesis. Journal of Bone Mineral Research. 2005;20(5):840–7.

    CAS  PubMed  Article  Google Scholar 

  63. 63.

    Saftig P, Hunziker E, Everts V, Jones S, Boyde A, Wehmeyer O, Suter A, von Figura K: Functions of cathepsin K in bone resorption. In: Cellular Peptidases in Immune Functions and Diseases 2. Springer; 2002: 293–303.

  64. 64.

    Levy-Apter E, Finkelshtein E, Vemulapalli V, Li SS-C, Bedford MT, Elson A. Adaptor protein GRB2 promotes Src tyrosine kinase activation and podosomal organization by protein-tyrosine phosphatase ϵ in osteoclasts. J Biol Chem. 2014;289(52):36048–58.

    CAS  PubMed  PubMed Central  Article  Google Scholar 

  65. 65.

    Paternoster L, Lorentzon M, Lehtimäki T, Eriksson J, Kähönen M, Raitakari O, Laaksonen M, Sievänen H, Viikari J, Lyytikäinen L-P. Genetic determinants of trabecular and cortical volumetric bone mineral densities and bone microstructure. PLoS Genet. 2013;9(2):e1003247.

    CAS  PubMed  PubMed Central  Article  Google Scholar 

  66. 66.

    Franz-Odendaal TA, Hall BK, Witten PE. Buried alive: how osteoblasts become osteocytes. Dev Dyn. 2006;235(1):176–90.

    CAS  PubMed  Article  Google Scholar 

  67. 67.

    Alford AI, Hankenson KD. Matricellular proteins: extracellular modulators of bone development, remodeling, and regeneration. Bone. 2006;38(6):749–57.

    CAS  PubMed  Article  Google Scholar 

  68. 68.

    Fratzl P, Gupta H, Paschalis E, Roschger P. Structure and mechanical quality of the collagen–mineral nano-composite in bone. J Mater Chem. 2004;14(14):2115–23.

    CAS  Article  Google Scholar 

  69. 69.

    Zhou X, Sun D, Guang X, Ma S, Fang X, Mariotti M, Nielsen R, Gladyshev VN, Yang G. Molecular footprints of aquatic adaptation including bone mass changes in cetaceans. Genome Biology and Evolution. 2018.

  70. 70.

    Volk SW, Shah SR, Cohen AJ, Wang Y, Brisson BK, Vogel LK, Hankenson KD, Adams SL. Type III collagen regulates osteoblastogenesis and the quantity of trabecular bone. Calcif Tissue Int. 2014;94(6):621–31.

    CAS  PubMed  PubMed Central  Article  Google Scholar 

  71. 71.

    Yen J-L, Lin S-P, Chen M-R, Niu D-M. Clinical features of Ehlers-Danlos syndrome. J Formos Med Assoc. 2006;105(6):475–80.

    PubMed  Article  Google Scholar 

  72. 72.

    Horton WA, Dwyer C, Goering R, Dean DC. Immunohistochemistry of types I and II collagen in undecalcified skeletal tissues. Journal of Histochemistry & Cytochemistry. 1983;31(3):417–25.

    CAS  Article  Google Scholar 

  73. 73.

    Termine JD, Kleinman HK, Whitson SW, Conn KM, McGarvey ML, Martin GR. Osteonectin, a bone-specific protein linking mineral to collagen. Cell. 1981;26(1):99–105.

    CAS  PubMed  Article  Google Scholar 

Download references


We are particularly grateful to Dr. Kaiya Zhou, Mr. Xinrong Xu and all of the students who have ever studied or are studying at NNU for their assistance in the skeleton sample collection. Many thanks to Meixiu Chen and Kui Li for helping to prepare the bone thin sections. We thank Aurore Canoville (University of North Carolina, Raleigh, North Carolina, U.S.), Shoji Hayashi (Okayama University of Science, Okayama, JAPAN), Michel Laurin and Alexandra Houssaye (Muséum national d’histoire naturelle, France), who have agreed or given bone images for the present study. We thank Xuanmin Guang and Wencai Jie for helping learning with Perl and R coding. We are especially grateful to Siming Ma and Lin Zhang for giving advice on the statistical analyses. Many thanks are also given to Yunxia Yang, Zhengfei Wang, Kangli Zhu, Ran Tian, Zepeng Zhang, Simin Chai, Yuan Mu and Yulin Gai for support and discussions.


This research was financially supported and designed by the key project of the National Natural Science Foundation of China (NSFC) (grant no. 31630071 to GY), the National Key Program of Research and Development, Ministry of Science and Technology (grant no. 2016YFC0503200 to GY and SX); The National Science Fund for Distinguished Young Scholars (grant no. 31325025 to GY) and the Priority Academic Program Development of Jiangsu Higher Education Institutions (PAPD) to GY and SX, and the Natural Science Foundation of Jiangsu Province of China (grant number BK20141449) to SX supported the data collection and analysis.

Author information




GY and SX designed the study. DS performed sample collection and data analysis, then prepared the draft of the manuscript. XZ participated in the design of the methods of PGLS and the writing of the manuscript. ZY helped analyze the data. IS put forward valuable opinions of and suggestions of result interpretation, assisted with manuscript editing, and polished the paper. All authors read and approved the final manuscript.

Corresponding authors

Correspondence to Inge Seim or Guang Yang.

Ethics declarations

Ethics approval and consent to participate

This study was carried out in full compliance with the laws of the People’s Republic of China. No specific permits were required for our field investigation. Only cadaver material was used in this study and the animals were not killed for this study.

Consent for publication

Not applicable.

Competing interests

The authors declare that they have 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

Additional file 1:

 Table S1-3. Datasets of rib, humerus and vertebrae compiled in the present study, respectively. Containing the list of species sampled, the collection number, the habitat, and the bone histological parameters measured for each species. Table S4. List of ‘root-to-tip’ dN/dS values of each species for 348 bone function-associated genes.

Additional file 2:

 Table S5. Data from phylogenetically informed ANOVA. Summary of pairwise comparison results of bone quantitative indices in different habitats using the phytools package in R. a. data from rib. b. data from humerus. Figure S1. Linear regressions result of MD and bone parameter (S, Cc, Cp) from rib; Linear regression result of CL and Vpc1 from vertebrae. Figure S2. A phylogeny of 24 marine mammals with corresponding bone section and the max diving depth for each species. The section images are Tursiops truncatus PL15-0145-28, Delphinus capensis PL15-0145-07, Stenella attenuata PL15-0145-14, Sousa chinensis PL15-0145-16, Globicephala melas AMNH 215271, Pseudorca crassidens PL15-0145-03, Lissodelphis borealis AMNH 31422, Orcinus orca AMNH 34261, Neophocaena asiaeorientalis PL15-0145-21, Phocoena phocoena MHNL 50001046, Delphinapterus leucas AMNH 34936, Monachus monachus MHNL 50001018, Inia goeffrensis AMNH 209101, Lipotes vexillifer PL15-0145-18, Mesoplodon densirostris AMNH 139931, Balaenoptera acutorostrata scammoni PL15-0145-08, Zalophus californianus AMNH 63946, Eumetopias jubatus AMNH 38400, Arctocephalus pusillus AMNH 81701, Odobenus rosmarus MHNL 50001014, Cystophora cristata AMNH 184659, Phoca vitulina MHNL 50001020, Monachus monachus MHNL 50001018, Enhydra lutris MHNL 50001023. Figure S3. Box plot of six habitat categories for five bone indices from ribs. Figure S4. The linear regressions of max/average diving depth with Cp and Cg for 24 aquatic mammals.

Additional file 3:

 Table S6. List of positively selected genes (PSGs) and divergently selected genes (DSGs) in cetaceans identified from 348 bone function-associated genes. Table S7. List of divergently selected genes (DSGs) in cetaceans identified from 348 bone function-associated genes.

Additional file 4:

 Table 8-10. Result of regression analyses between ‘root-to-tip’ dN/dS and bone histological parameters for rib, humerus and vertebrae, respectively.

Additional file 5:

 Table S11. Enrichment pathways for bone functional genes which ‘root-to-tip’ dN/dS significantly correlated with bone histological parameters. Table S12. Enrichment GO categories for bone functional genes which ‘root-to-tip’ dN/dS significantly correlated with bone histological parameters.

Rights and permissions

Open Access This article is distributed under the terms of the Creative Commons Attribution 4.0 International License (, which permits unrestricted use, distribution, and reproduction in any medium, provided you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons license, and indicate if changes were made. The Creative Commons Public Domain Dedication waiver ( applies to the data made available in this article, unless otherwise stated.

Reprints and Permissions

About this article

Verify currency and authenticity via CrossMark

Cite this article

Sun, D., Zhou, X., Yu, Z. et al. Accelerated evolution and diversifying selection drove the adaptation of cetacean bone microstructure. BMC Evol Biol 19, 194 (2019).

Download citation


  • Bone microstructure
  • Cetaceans
  • Adaptive evolution
  • Phylogenetic comparative analyses