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

Background 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. Results 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. Conclusions Our results suggest that evolutionary pressures have shaped the bone microstructure of cetaceans, to facilitate life in diverse aquatic environments.


Background
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 calciumphosphate 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 boneremodeling 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 labthe false killer whale (Pseudorca crassidens), long-beaked common dolphin (Delphinus capensis), minke whale (Balaenoptera acutorostrata scammony), pantropical spotted dolphin (Stenella 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 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].
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).
The 'root-to-tip' d N /d S , defined as the average value of accumulated d N /d S 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 d N /d S [34]. Regression between 'root-to-tip' d N /d S 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.
A comparison of 'root-to-tip' d N /d S 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' d N /d S and Cp and a positive correlation between 'root-to-tip' d N /d S and bone indices Cc and S, (Fig. 4b). Thus, genes associated with osteoclast function have evolved rapidly in aquatic mammals.
The association of 'root-to-tip' d N /d S 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 d N /d S 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 boneas 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' d N /d S 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' d N /d S 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' d N /d S of TNFRSF11B (osteoprotegerin; also known as OPG) positively correlated with Cc in both rib (P value.robust < 3.89 × 10 − 06 , r 2 = 0.91) and humerus (P value.robust < 2.10 × 10 − 06 , r 2 = 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 , r 2 = 0.75). In Fig. 3 Mapping of genes with 'root-to-tip' d N /d S significantly correlated with bone microstructures or under positive and divergent selections to the KEGG osteoclast differentiation pathway GWAS meta-analyses these genes were associated with volumetric bone mineral density, suggesting that the RANK-RANKL-OPG axis affects the skeletonat 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' d N /d S negatively correlated with Cc (P value.robust = 3.93 × 10 − 05 , r 2 = 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, 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' d N /d S ). 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' d N /d S ). 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 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.

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

Supplementary information
The online version of this article (https://doi.org/10.1186/s12862-019-1509-x) contains supplementary material, which is available to authorized users.
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' d N / d S 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.  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.