Odontogenic ameloblast-associated (ODAM) is inactivated in toothless/enamelless placental mammals and toothed whales

Abstract

Background

The gene for odontogenic ameloblast-associated (ODAM) is a member of the secretory calcium-binding phosphoprotein gene family. ODAM is primarily expressed in dental tissues including the enamel organ and the junctional epithelium, and may also have pleiotropic functions that are unrelated to teeth. Here, we leverage the power of natural selection to test competing hypotheses that ODAM is tooth-specific versus pleiotropic. Specifically, we compiled and screened complete protein-coding sequences, plus sequences for flanking intronic regions, for ODAM in 165 placental mammals to determine if this gene contains inactivating mutations in lineages that either lack teeth (baleen whales, pangolins, anteaters) or lack enamel on their teeth (aardvarks, sloths, armadillos), as would be expected if the only essential functions of ODAM are related to tooth development and the adhesion of the gingival junctional epithelium to the enamel tooth surface.

Results

We discovered inactivating mutations in all species of placental mammals that either lack teeth or lack enamel on their teeth. A surprising result is that ODAM is also inactivated in a few additional lineages including all toothed whales that were examined. We hypothesize that ODAM inactivation is related to the simplified outer enamel surface of toothed whales. An alternate hypothesis is that ODAM inactivation in toothed whales may be related to altered antimicrobial functions of the junctional epithelium in aquatic habitats. Selection analyses on ODAM sequences revealed that the composite dN/dS value for pseudogenic branches is close to 1.0 as expected for a neutrally evolving pseudogene. DN/dS values on transitional branches were used to estimate ODAM inactivation times. In the case of pangolins, ODAM was inactivated ~ 65 million years ago, which is older than the oldest pangolin fossil (Eomanis, 47 Ma) and suggests an even more ancient loss or simplification of teeth in this lineage.

Conclusion

Our results validate the hypothesis that the only essential functions of ODAM that are maintained by natural selection are related to tooth development and/or the maintenance of a healthy junctional epithelium that attaches to the enamel surface of teeth.

Background

The secretory calcium-binding phosphoprotein (SCPP) gene family contains 23 members in the human genome and achieved its present diversity through an extensive series of gene duplication events [1, 2]. In mammals, different SCPP genes are expressed in mineralized tissues, mammary glands, and salivary glands [2]. Among the genes that are expressed in mineralized tissues of vertebrates, some are important for bone and/or dentin formation whereas others are critical for the synthesis of enamel (tetrapods) and/or enameloid (bony fishes) [2]. The gene that encodes odontogenic ameloblast-associated (ODAM) is a member of the SCPP family and is thought to have originated as far back as the common ancestor of tetrapods and bony fishes more than 450 million years ago [2,3,4]. Early in tooth development, ODAM is highly expressed in the enamel organ that gives rise to ameloblasts that produce enamel. Later in tooth development, ODAM is expressed in the junctional epithelium, which is a specialized epithelium that attaches the soft tissue of the gingiva (gums) to the enamel surface of the tooth and protects against bacteria that promote periodontal disease [5,6,7,8,9]. In this capacity, ODAM plays a critical role in the first line of defense against bacterial invasion [10]. ODAM has also been shown to upregulate the expression of the MMP20 gene [11], which encodes a matrix metalloproteinase required for proper enamel development, and to promote the nucleation of hydroxyapatite crystals [12]. Young ODAM knockout mice show no abnormalities in enamel volume, density, and organization, but in older individuals the junctional epithelium shows significant detachment and apical downgrowth [13]. Beyond its expression in ameloblasts and the junctional epithelium, ODAM is expressed in many non-dental tissues in humans [8, 14]. Of these, expression levels are strongest in salivary glands, mammary glands, and the trachea [8]. ODAM is also expressed in carcinomas of the stomach, breast, lung, and ovary [14, 15].

Previous studies have shown that one or more of nine tooth-related genes (ACP4, AMBN, AMEL, AMTN, C4orf26 (= ODAPH), DSPP, ENAM, KLK4, MMP20) are inactivated in a variety of toothless vertebrates that have been investigated including birds [16,17,18,19], turtles [20], and several mammalian lineages comprising baleen whales, pangolins, anteaters, and Steller’s sea cow [19, 21,22,23,24,25,26,27,28]. Many of these genes are also inactivated in mammals with enamelless teeth such as pygmy and dwarf sperm whales, aardvarks, sloths, and armadillos [19, 23, 24, 28,29,30]. The widespread inactivation of these genes in edentulous and enamelless vertebrates suggests that their only essential functions are related to tooth development even if they are sometimes expressed in non-dental tissues [19].

Given its important role in amelogenesis, ODAM has emerged as yet another gene that should be screened in toothless and enamelless vertebrates for evidence of inactivating mutations. If ODAM is consistently inactivated in vertebrates that lack teeth or the outer enamel covering on their teeth, then ODAM’s only essential functions would appear to be tooth-related even though ODAM is expressed in other tissues [14]. By contrast, if ODAM gene sequences are intact in toothless and enamelless vertebrates this would provide evidence for pleiotropic, non-dental functions that affect fitness and are protected by natural selection.

Here, we test the hypothesis that the main function of ODAM is related to enamel formation by investigating whether this gene remains under purifying selection in placental mammals that have lost their teeth or the outer enamel covering on their teeth. To this end, we compiled and screened complete protein coding sequences of the ODAM gene for inactivating mutations in 165 placental mammal species comprising representatives of all extant placental orders. In support of our hypothesis, we found that all lineages without teeth (baleen whales, pangolins, anteaters) or without enamel on their teeth (aardvark, sloths, armadillos) have inactivated ODAM. Finally, selection analyses to estimate when purifying selection on ODAM was relaxed in lineages with inactivated ODAM provide new insights into when teeth or their outer enamel caps were lost in these lineages.

Methods

Taxon sampling

Taxon sampling included 165 placental mammal species, of which 135 are from assembled genomes in public databases, 13 are from unassembled sequence reads in the Sequence Read Archive (SRA) or European Nucleotide Archive (ENA), eight are from a combination of assembled genomes and unassembled sequence reads, and nine are from new Illumina HiSeq2500, HiSeq4000, HiSeqX10, or NovaSeq data (Additional file 1: Table S1). Taxon sampling included 12 taxa from Afrotheria, ten taxa from Xenarthra, 61 taxa from Euarchontoglires, and 82 taxa from Laurasiatheria. Among the afrotherians, three taxa (Mammut americanum [American mastodon], Mammuthus primigenius [woolly mammoth], Palaeoloxodon antiquus [straight-tusked elephant]) are extinct. Toothless taxa included two vermilinguans (anteaters), two pholidotans (pangolins), and seven mysticetes (baleen whales); taxa that lack enamel on their teeth included one tubulidentate (aardvark), two folivorans (sloths), and six cingulatans (armadillos).

BLAST searches and alignments

Genomic sequences encoding ODAM were downloaded from NCBI for representatives of Euarchontoglires (Homo sapiens), Laurasiatheria (Bos taurus), and Afrotheria (Trichechus manatus). These sequences were employed in BLAST searches against other placental mammals in NCBI’s 'RefSeq Genome' and 'Whole-genome shotgun contigs' databases. We used megablast for highly similar sequences and the more sensitive blastn for less similar sequences when megablast searches failed to return results. Complete coding sequences and intervening introns were imported into Geneious 11.1.5 [31] and aligned with MAFFT [32, 33] with minor adjustments by eye. For five xenarthrans (Cabassous unicinctus, Tolypeutes matacus, Euphractus sexcinctus, Choloepus didactylus, Tamandua tetradactyla), we imported Discovar de novo (https://software.broadinstitute.org/software/discovar/blog/) genome assemblies into Geneious [12] and performed local BLAST searches (discontiguous megablast) using reference xenarthran sequences (Dasypus novemcinctus, Choloepus hoffmanni). With a few exceptions (see below), exon-intron boundaries are conserved among placental mammals. In cases where there were suspected assembly gaps with missing exons, we used sequences from the closest available taxon (exon plus ~ 100 bp of each flanking intron) to perform additional BLAST searches with megablast against the Sequence Read Archive (SRA). Sequences from SRA BLAST searches were imported into Geneious and mapped to the closest reference sequence to assemble missing exons and their flanking intron sequences. A similar strategy was employed to retrieve complete protein-coding sequences for ODAM in taxa with unassembled genomes in the SRA. We also used a mapping to reference approach in Geneious to retrieve ODAM sequences from unpublished Illumina data for Choeropsis liberiensis (pygmy hippopotamus) and three xenarthrans. For C. liberiensis, we mapped reads to the ODAM sequence for Hippopotamus amphibius (river hippopotamus). For the three xenarthrans (Chlamyphorus truncatus, Calyptophractus retusus, Cyclopes didactylus), we imported short reads into Geneious and mapped them to various xenarthran reference sequences. Individual exon sequences (and flanking introns) from each SRA taxon were concatenated into a single sequence with Geneious and aligned to the global alignment.

Inactivating mutations

We inspected the final alignment (165 taxa) for inactivating mutations including exon deletions, frameshift mutations, altered start and stop codons, premature stop codons, and splice site mutations. We also used the alignment program MACSE v2 [34, 35] to crosscheck the results of manual inspections for frameshifts and premature stop codons in coding sequences. Parsimony optimizations with delayed transformation (deltran) were performed with PAUP* 4.0a150 [36] to map inactivating mutations to branches of the phylogeny of placental mammals.

Phylogenetic analyses

RAxML 8.2.10 [37, 38] on CIPRES [39] was used to estimate a maximum likelihood tree for the complete protein-coding sequence alignment. Rapid bootstrap analysis (500 pseudoreplications) [40] and a search for the best tree were performed in a single run. These analyses were performed with a GTR + Γ model of sequence evolution.

Selection analyses

We first reconstructed the evolution of the dN/dS ratio of ODAM across the placental mammal phylogeny using the Bayesian approach implemented in Coevol 1.4b [41] on the complete alignment of 165 sequences containing 279 aligned codons. We used the dsom estimation procedure that jointly estimates branch specific dN/dS ratios, divergence times, body sizes, generation times, and ages at sexual maturity modeled as a multivariate Brownian diffusion process [42]. DN/dS analyses with Coevol were performed to detect and visualize patterns of dN/dS variation across a tree. The assumed topology was the ODAM gene tree after making several rearrangements to correct for discrepancies between the gene tree and species tree [43,44,45,46] as follows: Chaetophractus sister to other Chlamyphoridae; Orycteropus sister to Afrosoricida + Macroscelidea; Solenodon sister to other Eulipotyphla; Sus sister to Ruminantia + Whippomorpha; Physeter sister to other Odontoceti; Pholidota sister to Carnivora; Rhinolophus + Hipposideros sister to Megaderma; and Nomascus sister to Hominidae. We assumed fossil calibrations from previous studies [44, 46, 47] and extracted life history traits from PanTHERIA [48]. We set the prior on the root node to 97 Ma with a standard deviation of 20 Ma following the dating of Meredith et al. [44]. We ran two independent MCMC for a total of 5000 cycles sampling parameters every cycle. After checking for convergence by monitoring the effective sample size of the different parameters using the tracecomp command, we excluded the first 500 points of each MCMC as burnin, and made inferences from the remaining 4500 sampled points of each chain.

Additional dN/dS analyses were performed with the codeml program of PAML 4.4 [49]. Codeml does not employ the autocorrelative model of Brownian motion [50] that is incorporated into Coevol. For computational efficiency and tractability, codeml analyses were performed with 96 taxa after pruning all but one representative for mammalian families with ODAM sequences that are intact for all exemplars included in our study. We used 15 dN/dS categories based on reconstructions of ancestral states for tooth/enamel loss [23, 24] and the occurrence of frameshift mutations, premature stop codons, and/or exon deletions in independent lineages. We ignored taxa (Hyaena hyaena, Elaphurus davidianus, Eulemur flavifrons [see below]) that contained just a single splice site mutation except for Physeter macrocephalus. This taxon is of interest because the next possible donor splice site (GT) for intron 4 is 75 bp downstream from the canonical splice site and has a very low MaxEntScan score of − 17 (http://genes.mit.edu/burgelab/maxent/Xmaxentscan_scoreseq.html [51]). (For reference, splice sites are predicted to be functional when they have MaxEntScan scores of > 3 [52]). Moreover, all other toothed whales have one or more frameshift mutations and/or premature stop codons in ODAM. The 15 categories that we employed in dN/dS analyses are as follows: one category for branches leading to intact sequences that lack inactivating mutations; one category for pseudogenic branches that postdate tooth or enamel loss (parsimony deltran optimization) and/or the occurrence of an inactivating mutation on an earlier branch; nine categories for transitional branches that record the first inactivating mutation in ODAM (Dasypus novemcinctus; stem Chlamyphoridae; Orycteropus afer; stem Pholidota; Neomonachus schauinslandi; Lipotes vexillifer; stem to Delphinidae + Monodontidae + Phocoenidae; Physeter macrocephalus; Callithrix jacchus); one category for the stem Pilosa branch where enamel was lost [23]; one category for the stem Mysticeti branch where teeth were lost [24]; one category for the stem branch leading to proboscideans with missing exon 6; and one category for proboscidean branches that post-date the loss of exon 6. We employed the same topological corrections for the species tree as detailed above for the Coevol analysis. Codeml analyses were performed with two different codon frequency models, CF1 and CF2. Codon frequencies are estimated from mean nucleotide frequencies across all three codon positions in CF1 and from mean nucleotide frequencies at each of the individual codon positions in CF2. All frameshift insertions were removed from the alignment prior to codeml analyses, as were in-frame insertions that are unique to one or a few taxa. Also, premature stop codons were recoded as missing data. The final alignment for codeml analyses contained 837 aligned positions (279 codons).

Estimation of ODAM inactivation times

We used equations from Meredith et al. [23] to estimate when ODAM was inactivated in different placental lineages. Calculations were performed with dN/dS values that were obtained using two different codon models in codeml (CF1, CF2), fixed (1.0) versus estimated values for the dN/dS value on fully pseudogenic branches, and equations that allow for one versus two synonymous substitution rates [23]. Divergence times for relevant nodes in these calculations were taken from Springer et al. [53] for Callithrix to Aotus, Delsuc et al. [45] for all divergences in Xenarthra, Gaubert et al. [54] for Manis pentadactyla to M. javanica, and Foley et al. [46] for all other divergence dates.

Screen for convergently inactivated genes in edentulous and enamelless mammals

We performed a genomic screen for convergently inactivated genes based on a Forward Genomics approach with a multispecies genome alignment [28]. Briefly, this screen is based on the percentage of the reading frame of the protein-coding gene that is intact and is calculated based on the relative position of inactivating mutations and the partial and/or complete loss of entire exons. This percentage is hereafter denoted as %intact. We considered 61 placental mammals for which inactivating mutations have been computed previously [28]. These species include the edentulous Balaenoptera acutorostrata scammoni and Manis pentadactyla and the enamelless Orycteropus afer and Dasypus novemcinctus, representing four independent lineages. In contrast to a previous screen [28], we performed a more sensitive search by classifying a gene as inactivated if %intact is < 80% and classifying a gene as intact if %intact is ≥90%. We then screened for genes that are inactivated in at least two of the four edentulous or enamelless species, but are intact in at least 90% of species that have tooth enamel. The resulting candidate list of the screen was then sorted (largest to smallest) by the number of species without teeth or enamel that exhibit an inactivated gene (%intact < 80%) and, in case of a tie, sorted (smallest to largest) by the number of species with tooth enamel that do not have an intact version of the gene (%intact < 90%). Finally, we used a phylogenetic generalized least squares approach [55] to compute a phylogeny-corrected P-value for the association between %intact and the presence/absence of enamel.

Results

Multiple alignment and phylogenetic analysis

The 165-taxon alignment for the complete protein-coding sequence of ODAM is 1177 nucleotides including frameshift insertions and is available at TreeBASE (http://purl.org/phylo/treebase/phylows/study/TB2:S23531). The shape parameter of the gamma distribution was estimated at 1.98 by RAxML. Figure 1 shows the maximum likelihood gene tree (ln L = − 27,926.663886). This gene tree was not rooted with a marsupial, but is consistent with the monophyly of the four major clades of placental mammals (Xenarthra, Afrotheria, Laurasiatheria, Euarchontoglires) [56,57,58,59,60,61,62,63,64]. Other superordinal groups that were recovered on the ODAM gene tree include Afroinsectiphilia, Paenungulata, Variamana, Euungulata, Glires, and Euarchonta. Thirteen of 19 placental orders included more than one species and in all cases these orders were recovered as monophyletic. Several features of the ODAM gene tree are in minor conflict with generally accepted species trees including the placement of Physeter as the sister group to mysticetes rather than other odontocetes [65, 66], an association of Pholidota with Chiroptera rather than Carnivora [44], a sister group relationship between Hipposideridae + Rhinolophidae and yangochiropteran bats instead of with other rhinolophoid bats [43, 67], an association of chlamyphorine armadillos (Chlamyphorus, Calyptophractus) with Chaetophractus rather than with tolypeutines (Tolypeutes, Cabassous) [45], a basal split between Condylura and other eulipotyphlans instead of between Solenodon and other eulipotyphlans [44, 68, 69], and the nonmonophyly of Elephantidae with the mammutid Mammut americanum nested inside of this clade instead of sister to this clade [70, 71].

Fig. 1

Maximum likelihood phylogram (ln L = − 27,926.663886) based on ODAM protein-coding sequences for 165 placental mammal species. Pink boxes indicate taxa with inactivating mutations in ODAM; green boxes indicate proboscideans that are missing exon 6 of ODAM. The tree is rooted between Atlantogenata and Boreoeutheria. Some clades for which all of the constituent taxa have intact ODAM sequences were collapsed, and the length of the Glires clade branch was halved to improve aesthetics (see TreeBASE [http://purl.org/phylo/treebase/phylows/study/TB2:S23531] for ML phylogram with all 165 taxa)

Inactivating mutations

One or more inactivating mutations (frameshifts, altered start and stop codons, premature stop codons, splice site mutations) were discovered in all placental taxa without teeth or without enamel on their teeth. Table 1 provides a complete list of inactivating mutations in these taxa and examples of inactivating mutations are shown in Fig. 2. Shared inactivating mutations were discovered in Vermilingua, Choloepus, Chlamyphoridae, Tolypeutinae, Manis, Balaenidae, Balaenopteroidea, and in the two minke whales (Balaenoptera acutorostrata and B. bonaerensis). Among edentulous forms, taxa with the largest number of inactivating mutations are Manis javanica (15), Manis pentadactyla (14), and Tamandua tetradactyla (14), whereas the fewest inactivating mutations are found in Balaena mysticetus (1) and Balaenoptera musculus (1). Among taxa with enamelless teeth, Orycteropus afer and Tolypeutes matacus each have four inactivating mutations whereas Choloepus didactylus only has one inactivating mutation. However, it should be noted that only 8 of 10 exons had BLAST results for C. didactylus. An association between ODAM inactivation and the loss of teeth or enamel is further supported by a genomic screen for genes that have inactivating mutations and an abolished reading frame in edentulous and enamelless taxa, which uncovered ODAM at rank 7 together with other tooth-related genes (Additional file 2: Table S2).

Table 1 Inactivating mutations and whole exon deletions in the ODAM gene of toothless, enamelless, and other placental mammals

Fig. 2

Examples of inactivating mutations in ODAM. a) Orycteropus afer (aardvark). b) Xenarthra (anteaters, sloths, armadillos). c) Manis spp. (pangolins). d) Cetacea (whales, dolphins, porpoises). Green boxes denote frameshift insertions; blue boxes denote frameshift deletions; and purple boxes denote premature stop codons. Loxodonta, Procavia, Hyaena, Canis, and Hippopotamus have intact ODAM sequences and were included in different alignment panels to provide context for inactivating mutations

In addition to inactivating mutations in edentulous and enamelless taxa, we also discovered inactivating mutations in all eight representatives of Odontoceti including inactivating mutations that are shared by members of the clade Delphinoidea (Phocoenidae [3 spp.], Monodontidae [1 sp.], and Delphinidae [2 spp.]). The number of inactivating mutations among different odontocete species ranges from one in Physeter macrocephalus (splice site mutation) to four in the three phocoenid species (Phocoena phocoena, Neophocaena phocaenoides, N. asiaeorientalis).

Beyond these odontocetes, additional taxa with enamel-capped teeth had putative inactivating mutations. The phocid (seal) Neomonachus schauinslandi has one 2-bp frameshift deletion, one premature stop codon, an acceptor splice site mutation (intron 2 acceptor, AG = > AA), and a putative deletion of exon 5. However, genomic Illumina data are not available for Neomonachus, and we were unable to validate these mutations. By contrast, the New World primate Callithrix jacchus has a premature stop codon in exon 5 and an acceptor splice site mutation (AC) in intron 4, both of which are confirmed by SRA data. The mustelid Enhydra lutris has a 1-bp frameshift insertion in exon 1. Here, Illumina data for three individuals of Enhydra (SRR5768046, SRR5768052, SRR6450476) suggest that this deletion either has allelic variation or that there is an additional paralog of exon 1 in the genome of this species. Finally, there are splice site mutations (relative to canonical splice sites) in the strepsirrhine primate Eulemur flavifrons (intron 7 acceptor, AG = > TT), in the carnivoran Hyaena hyaena (intron 8 donor, GT = > TT), and in the deer Elaphurus davidianus (intron 5 donor, GT = > CT). However, all of these splice site mutations can be accommodated by alternative splice sites that result in slightly shorter exons (3 bp shorter in Elaphurus, 6 bp shorter in Hyaena, 18 bp shorter in Eulemur), and such evolutionary splice site shifts provide no indication for gene loss [72].

Finally, there were no BLAST results for exon 6 of ODAM in three elephantids (Mammuthus primigenius, Palaeoloxodon antiquus, Loxodonta cyclotis) and the mammutid Mammut americanum (Fig. 3). By contrast, this exon is present in the elephantids Loxodonta africana and Elephas maximus.

Fig. 3

Map to reference coverage of Illumina reads onto a contiguous block of Loxodonta africana (African savannah elephant) ODAM that begins with intron 5 and ends with exon 10. Map to reference results suggest that exon 6 is retained in Elephas maximus (Asian elephant) but is missing in several other proboscideans (Loxodonta cyclotis [African forest elephant], Mammuthus primigenius [woolly mammoth], Palaeoloxodon antiquus [straight-tusked elephant], Mammut americanum [American mastodon])

Selection analyses

The joint Bayesian reconstruction of dN/dS across the placental phylogeny highlights large variation of selection pressure among lineages. Mean dN/dS values for functional branches range from 0.36 in the kangaroo rat Dipodomys ordii to 0.70 in the bonobo Pan paniscus and are compatible with purifying selection (Fig. 4). At the ordinal level, rodents and bats show the lowest average dN/dS values of 0.41 and 0.48, respectively, whereas functional cetartiodactyl and primate sequences exhibit larger mean values of 0.59 and 0.60, respectively. These results illustrate the influence of different life-history traits. A consequence of relaxed selection constraints is that the non-functional sequences show elevated dN/dS values for pilosans (0.79), cingulatans (0.81), aardvark (0.73), pangolins (0.70), odontocetes (0.75) and mysticetes (0.87). However, the Hawaiian monk seal (Neomonachus schauinslandi) and the common marmoset (Callithrix jacchus), which present inactivating mutations in their ODAM sequences, show only moderately elevated dN/dS values of 0.58 and 0.54, respectively. Transitional branches where enamel loss was inferred based on ancestral reconstructions [23, 24] all have elevated dN/dS values: stem Pilosa (0.75), stem Chlamyphoridae (0.81), stem Pholidota (0.70), and stem Mysticeti (0.82).

Fig. 4

Bayesian reconstruction of dN/dS for 165 ODAM sequences across the placental phylogeny. The variation of dN/dS was jointly reconstructed with divergence times while controlling the effect of three life-history traits (body mass, longevity, and maturity). Asterisks at the tips of terminal branches indicate non-functional sequences (pseudogenes). Stars indicate branches on which shared inactivating mutations were inferred in toothless or enamelless clades. The tree is rooted between Atlantogenata and Boreoeutheria. Placental orders are delimited to the right of species tree tips

The results of selection (dN/dS) analyses using codeml are summarized in Table 2. The dN/dS values for functional branches that lead to taxa with intact protein-coding sequences for ODAM and enamel-capped teeth are ~ 0.49 and ~ 0.52 with codon frequency model 1 (CF1) and codon frequency model 2 (CF2), respectively, confirming that ODAM has evolved overall under purifying selection. By contrast, the values for pseudogenic branches that post-date the first occurrence of an inactivating mutation(s) on an earlier branch are ~ 1.04 (CF1) and 1.11 (CF2). These values are close to the expected value of 1.00 for neutrally evolving pseudogene sequences. Transitional branches where enamel was lost based on ancestral reconstructions [23, 24] include stem Pilosa, stem Pholidota, stem Mysticeti, and Orycteropus. DN/dS values for these branches are ~ 0.12 (CF1) to ~ 0.14 (CF2) for stem Pilosa, ~ 0.93 (CF1) to ~ 1.0 (CF2) for stem Pholidota, and ~ 0.88 (CF1) to ~ 0.93 (CF2) for Orycteropus. The codeml value for Mysticeti is based on only two nucleotide substitutions, both of which are nonsynonymous, leading to parameter unidentifiability. Other transitional branches for enamelless species based on the occurrence of the first inactivating mutation(s) include stem Chlamyphoridae and Dasypus. The dN/dS values for stem Chlamyphoridae are ~ 1.01 (CF1) and ~ 1.10 (CF2), whereas those for Dasypus are ~ 0.96 (CF1) and ~ 1.02 (CF2).

Table 2 Results of selection (dN/dS=ω) analyses with codon frequency models 1 (CF1) and 2 (CF2) and 15 different branch categories

Transitional branches that exhibit the first inactivating mutation(s) among odontocetes all have high dN/dS values including Physeter (~ 0.90 [CF1], ~ 1.01 [CF2]), and Lipotes (~ 1.11 [CF1], ~ 1.17 [CF2]). The dN/dS value for the stem branch leading to Delphinidae + Phocoenidae + Monodontidae is based on 7.5–7.7 substitutions, all of which are nonsynonymous, and is undefined. Other transitional branches with putative inactivating mutations are Callithrix, Neomonachus, and the stem branch leading to three elephantids + Mammut americanum. Callithrix has an elevated dN/dS value (~ 0.85 [CF1], ~ 0.90 [CF2]). Neomonachus has dN/dS values of ~ 0.40 (CF1) to ~ 0.41 (CF2) that are below the background ratios for codeml. Finally, there are no inferred substitutions on the stem branch leading to proboscideans that are missing exon 6 so the codeml values for this branch can be ignored. The dN/dS values for crown proboscidean branches that are missing exon 6 are ~ 0.21 (CF1) and ~ 0.23 (CF2) and are below the background ratios of ~ 0.49 and ~ 0.52, respectively.

Inactivation times

Estimates of inactivation times for ODAM based on dN/dS ratios and equations in Meredith et al. [23] are provided in Table 3 for several transitional branches. The mean estimate for the inactivation of ODAM on the stem Pholidota branch is 64.15 Ma (range = 73.04–57.34 Ma) (Fig. 5). The mean inactivation estimate for Orycteropus ODAM, in turn, is 55.63 Ma (range 58.9–48.2 Ma). Among cingulatans, the mean inactivation time in Dasypus is 40.43 Ma (range 45.45–36.38 Ma) and is similar to the mean estimate of 45.29 Ma for the stem chlamyphorid branch (range from 45.45–44.73 Ma). Inactivation dates for Lipotes and the stem Delphinidae + Phocoenidae + Monodontidae branch were not calculated because dN/dS values for both branches always exceeded the estimates for pseudogenic branches. However, the elevated dN/dS values for these branches suggest that selection on ODAM was relaxed near the base of each branch or even earlier (see Discussion). The mean inactivation date for Callithrix ODAM is 13.12 Ma (range = 14.47–11.17 Ma).

Table 3 Estimated inactivation times, in millions of years, for select branches where the ODAM gene accumulated its first inactivating mutation

Fig. 5

Estimated inactivation times in ENAM versus ODAM. a) Manis spp. (pangolins). b) Orycteropus afer (aardvark). c) Pilosa (anteaters, sloths) and Dasypus novemcinctus (nine-banded armadillo). Inactivation dates are mean values based on eight different combinations of two different codon models in codeml (CF1, CF2), fixed (1.0) versus estimated values for the dN/dS value on fully pseudogenic branches, and equations that permit for one versus two synonymous substitution rates [23]. Manis tricuspis is missing from the ODAM data set and M. javanica was not included in Meredith et al.’s [23] ENAM data set

Discussion

Patterns of ODAM inactivation in placental mammals

Previous studies have documented inactivating mutations in nine different tooth-related genes in toothless and enamelless mammals (ACP4, AMBN, AMEL, AMTN, C4orf26, DSPP, ENAM, KLK4, MMP20). ODAM can now be added to this list based on inactivating mutations in both edentulous (Vermilingua, Pholidota, Mysticeti) and enamelless (Orycteropus, Folivora, Cingulata) clades. Indeed, inactivating mutations were discovered in all 19 edentulous and enamelless species that were investigated. Given the broad phylogenetic spread of mammalian lineages with inactivating mutations, these results suggest that the only essential functions of ODAM that are maintained by natural selection are related to tooth development and/or maintenance of the adhesion of the junctional epithelium to the tooth surface, even though ODAM expression has been reported in other tissues such as salivary gland, trachea, mammary gland, and lacrimal gland [14, 15]. The expression of ODAM in mammary and salivary glands may be explained by the location of this gene in the same gene expression neighborhood [73] as other SCPP genes that are expressed in mammary glands (CSN2, CSN3) and salivary glands (STATH).

Patterns of ODAM inactivation are generally consistent with ancestral reconstructions of tooth and enamel loss and previous studies of other tooth-related genes [19, 23, 24, 28, 30]. Our estimates for ODAM inactivation in Pholidota (73.04–57.34 Ma) are slightly older than Meredith et al.’s [23] estimates for ENAM inactivation (59.4–54.9) in this lineage. We recalculated inactivation dates for pangolin ENAM using the same divergence times and codon frequency models that were used for ODAM and the mean inactivation date is 57.7 Ma (Fig. 5). Estimates for ODAM and ENAM inactivation are both older than the oldest fossil pangolin, Eomanis waldi, which is ~ 47 Ma [23, 74], and suggest that even older fossil pholidotans that lack teeth or at least tooth enamel may be discovered.

Inactivation dates for ODAM are also older than inactivation dates for ENAM in the aardvark Orycteropus afer. Meredith et al. [23] estimated an inactivation date of 35.3–28.8 Ma for O. afer ENAM, and we obtained a similar estimate (~ 30.7 Ma) when we used the same divergence dates and codon frequency models as for ODAM (Fig. 5).

In Xenarthra, ancestral reconstructions of the presence or absence of enamel suggest that enamel was lost on the stem Pilosa branch [23]. This prediction is validated by two inactivating mutations in ENAM that are shared by Vermilingua (anteaters) and Folivora (sloths). In ODAM, we discovered inactivating mutations that are shared by the two vermilinguans included in our study, as well as a mutation that is shared by both folivorans, but did not find any mutations that are shared by all pilosans. The dN/dS value for the stem Pilosa branch (~ 0.12–0.14) is below the background dN/dS ratio of ~ 0.49–0.52 and suggests that ODAM was maintained by natural selection up until the split between Vermilingua and Folivora. We also re-ran codeml analyses with an added category for the stem Vermilingua and stem Folivora branches, which are the immediate descendant branches of stem Pilosa. The resulting dN/dS values (~ 1.12 [CF1], ~ 1.22 [CF2]) are slightly above the dN/dS values for the pseudogenic branch category (~ 1.04 [CF1], ~ 1.11 [CF2]) and suggest that ODAM has evolved neutrally on these branches. Together, the dN/dS values for stem Pilosa, stem Vermilingua, and stem Folivora suggest that selection on ODAM was relaxed very near the most recent common ancestor of Pilosa. By contrast, dN/dS analyses for ENAM suggest that selection was relaxed near the base of the stem pilosan branch at ~ 64.8 Ma (Fig. 5).

Also in Xenarthra, our estimates of inactivation times in Cingulata (armadillos) suggest that selection on ODAM was relaxed very soon after the split between Dasypodidae and Chlamyphoridae at 45.5 Ma [45]. Specifically, selection on ODAM was relaxed ~ 40 Ma in Dasypus and ~ 45 Ma in stem Chlamyphoridae. Meredith et al. [23] reported inactivating mutations in ENAM for four cingulatan genera that they investigated, but none that were shared by all chlamyphorids. In the case of ENAM, we estimated an inactivation date of ~ 45.5 Ma for this gene in Dasypus. By contrast with pilosans, where inactivating mutations in ENAM preceded inactivating mutations in ODAM, the opposite pattern occurs in chlamyphorids with an inactivating mutation in ODAM preceding the first inactivating mutations in ENAM (or at least the region of ENAM that was targeted by Meredith et al. [23]).

Whereas ENAM is required for enamel formation, ODAM inactivation in mouse only causes abnormalities in the junctional epithelium [13]. Thus, the more ancient loss of ODAM compared with ENAM in both Pholidota and Orycteropus afer suggests that selection on maintenance of the junctional epithelium was relaxed before enamel loss likely happened. The inactivation of ODAM but not ENAM in toothed whales that retain enamel also supports this hypothesis. However, Pilosa and Dasypus novemcinctus both show the opposite pattern where inactivation dates for ENAM are older than for ODAM (Fig. 5).

Beyond its inactivation in toothless and enamelless forms, we also discovered inactivating mutations or exon deletions in taxa with enamel-capped teeth. Two of these (odontocetes, some proboscideans) are discussed in separate sections below. Among the remaining taxa with mutations in ODAM, the phocid Neomonachus has four separate mutations (2 bp frameshift deletion, premature stop codon, splice site mutation, exon 5 deletion) and is the best candidate for inactivation of this gene among taxa with enamel-capped teeth. However, short read data for Neomonachus are not available in the SRA and we were unable to validate these mutations. Another phocid, Leptonychotes weddellii, also has problems with the annotated gene sequence, and like Neomonachus is missing exon 5 (as well as exons 6–8). SRA data are available for Leptonychotes and we were able to reconstruct all of the missing exons for this taxon. It remains unclear if the missing exon and other problems with Neomonachus ODAM are real or instead are assembly or annotation errors. Enhydra lutris has a 1 bp frameshift in exon 1 of the annotated sequence, but three different individuals with SRA data have both intact and frameshifted versions of this exonic region, which leaves open the possibility that the mutated form of exon 1 represents a second allele, or more likely, a paralogous locus. Callithrix jacchus (common marmoset) has an acceptor splice site mutation in intron 4, a premature stop codon in exon 5, and an elevated dN/dS ratio (~ 0.85–0.90). Curiously, marmosets exhibit a loss of lingual enamel and a hypertrophy of buccal enamel in their mandibular incisors [75]. However, similar phenotypes are found in the aye-aye (Daubentonia madagascariensis) and rodents, yet they uniformly retain an intact ODAM.

ODAM inactivation in odontocetes

It is perhaps surprising that ODAM contains inactivating mutations in all of the odontocetes that we investigated given that all of these taxa retain teeth with enamel. We hypothesize that ODAM inactivation in toothed whales may be related to the simplified outer enamel of their teeth. Living odontocetes are characterized by degenerative enamel that is thinner and less complex than the enamel of their archaeocete (stem cetacean) relatives [76,77,78]. In archaeocetes that have been investigated, there is an inner layer of enamel that is organized into Hunter-Schreger bands (HSB), which are decussating layers of prisms that increase the strength of enamel [79], and an outer layer of radial enamel. This pattern is common in large-bodied, terrestrial mammals and is thought to be a biomechanical adaptation for food processing and crack resistance [78]. By contrast, living odontocetes often swallow prey whole, without mastication, and have variable enamel structure with loss of the HSB in the inner layer. Among Delphinida, which includes all of the odontocetes in our study excepting the sperm whale (Physeter macrocephalus), the simplified enamel is typically comprised of an inner layer of radial enamel and an outer layer of aprismatic enamel [77, 78], although in some cases the entire enamel layer is aprismatic [76]. The enamel in Physeter has been characterized as being comprised of pseudoprisms. The simplified enamel microstructure of living odontocetes is associated with the transition from heterodont to homodont dentition wherein the upper and lower teeth are used to grasp and secure prey but not for mastication [77, 78]. Given that ODAM is expressed in the junctional epithelium, where the gingiva is in contact with the outer enamel surface of the tooth [5], perhaps ODAM inactivation in odontocetes is related to the simplified outer enamel surface if this surface is not well suited for ODAM-assisted adhesion. Indeed, ODAM is also inactivated in placental mammals that lack enamel on their teeth and instead present dentin on the outer surface (aardvarks, sloths, armadillos).

An alternate hypothesis is that ODAM inactivation in toothed whales is related to antimicrobial functions of the junctional epithelium that are altered in aquatic habitats. Specifically, the gingival junctional epithelium adheres tightly to the enamel surface of the tooth in terrestrial mammals and presents a first line of defense against invading bacteria [5]. However, the microbiomes of the oral cavity are expected to differ substantially in terrestrial versus aquatic environments, and toothed whales may require a different line of defense against microbes than is required by terrestrial mammals. The coding sequence of ODAM is intact in Trichechus manatus (West Indian manatee), which is the other fully aquatic mammal with teeth that was included in our analysis, but codeml analyses revealed an elevated dN/dS value on the Trichechus branch (CF1 = 1.60, CF2 = 1.89) when this branch was allowed to have its own category. Although these values are > 1 and suggestive of positive selection, there is no significant difference between a free ratio for Trichechus versus a dN/dS value of 1.0 for this branch (p = 0.33 with CF1, p = 0.18 with CF2). Thus, ODAM evolution on the Trichechus branch is also consistent with neutral evolution, as might be expected if a gene is evolving neutrally even though the first inactivating mutation has not yet been fixed. It is worth noting that manatees have a unique system of tooth replacement, hind molar progression, in which the molars march forward until they are worn down and replaced by new molars that emerge at the posterior end of the tooth row [80].

ODAM deletion and functionality in Proboscidea?

Our taxon sampling included six proboscideans, of which five (Loxodonta africana, L. cyclotis, Elephas maximus, Mammuthus primigenius, Palaeoloxodon antiquus) belong to Elephantidae and one belongs to Mammutidae (Mammut americanum). Exon 6 is putatively missing from four of these taxa based on the absence of sequencing reads that map to this region (Fig. 3). The assembled genome sequence for L. africana is also missing exon 6, but we were able to assemble this exon and its flanking intronic regions from SRA data where the Sanger sequencing-based assembly contains a string of Ns (assembly gap). The presence of exon 6 in L. africana and E. maximus, but not in other proboscideans, is unexpected given the phylogenetic relationships of these taxa. Mammutidae and Elephantidae diverged from each other ~ 30–20 Ma [70, 71], whereas the five elephantids share a much more recent common ancestor ~ 10–5 Ma [71, 81]. Figure 6 shows a time tree for the six proboscideans in our study along with presence/absence data for exon 6 of ODAM. Three possible explanations for the presence/absence of exon 6 in these taxa include incomplete lineage sorting (ILS), hybridization, and/or independent deletions of exon 6. Figure 6A shows how ILS can account for the presence/absence of exon 6 in proboscideans, whereas Fig. 6B illustrates how a combination of convergent loss and hybridization can explain the same pattern. Hybridization may play a role in the distribution of ODAM exon 6 in Elephantidae given recent evidence for extensive introgression among elephantids [71], but there is no evidence for hybridization between crown elephantids and Mammut americanum and this explanation seems unlikely. Instead, ILS or the independent deletion of exon 6 in Mammut americanum versus elephantids that lack this exon seems more probable. Is also remains to be determined if individual elephantid species will reveal presence/absence variation for exon 6 when more individuals are sampled, which would not be surprising in view of extensive hybridization among elephantids.

Fig. 6

Different hypotheses for the loss of ODAM exon 6 in some but not all proboscideans. a) Loss of ODAM exon 6 based on an ancestral polymorphism in the ancestor of Elephantidae and Mammutidae followed by incomplete lineage sorting of the two allelic variants in Recent and extinct proboscideans. b) One possible scenario for the loss of exon 6 of ODAM based on convergent loss in Mammut americanum and Palaeoloxodon antiquus followed by introgression of the sans exon 6 allele from P. antiquus to Loxodonta cyclotis and Mammuthus primigenius. Other scenarios are also possible based on directions of introgression within Elephantidae that were identified by Palkopoulou et al. [71]

It is also worth noting that the deletion of exon 6 seemingly has an ancient history that may trace as far back as the common ancestor of Elephantidae and Mammutidae if ILS is responsible for the presence/absence of this exon in different proboscideans. Exon 6 is 111 bp in elephantids that have this exon, so complete deletion of this exon will not interrupt the reading frame of the downstream exons (7–10). This leaves open the possibility of an intact version of ODAM based on nine of ten exons. The possibility of an altered yet functional version of ODAM is also supported by the absence of any frameshift mutations, premature stop codons, or splice site mutations in the remaining nine exons (1–5, 7–10) of the ODAM gene of proboscideans that are missing exon 6. Finally, all of the proboscidean ODAM sequences that lack exon 6 cluster together on the ODAM gene tree and the dN/dS value for this clade is only ~ 0.21–0.23, which is below the background value (~ 0.49–0.52) for branches leading to taxa with taxa with enamel-capped teeth and intact ODAM sequences for all ten exons. Thus, beyond the absence of exon 6, there is no evidence for relaxed selection in the ODAM gene of proboscideans that lack this exon. Unfortunately a 3D structure for ODAM is not available and it remains unclear if exon 6 encodes a region of the protein that is on the surface or embedded in the interior of this molecule.

Coevol versus codeml selection analyses

DN/dS analyses with Coevol and codeml were performed with different taxon sets and for different purposes. Coevol was used to implement a joint analysis that included all 165 taxa. This analysis was performed to visualize patterns of dN/dS variation across the placental tree and to examine the effects of life history traits on dN/dS values. Importantly, Coevol employs a form of rate smoothing through the incorporation of a Brownian motion model of continuous trait evolution. By contrast, codeml analyses were performed with 96 taxa and 15 different dN/dS categories that were determined a priori based on ancestral reconstructions of tooth/enamel loss [23, 24] and patterns of ODAM inactivation. Coevol and codeml results are generally in agreement with each other with lower dN/dS values for functional branches and higher values for transitional and fully pseudogenic branches. At the same time, dN/dS values for both transitional and pseudogenic branches are generally higher for codeml than for Coevol. For example, Coevol dN/dS values are < 1 for pseudogenic branches in clades that post-date the loss of teeth or enamel. By contrast, dN/dS values for the fully pseudogenic branch category are > 1 in codeml analyses (Table 2). Similarly, the dN/dS value for the transitional branch leading to Callithrix is higher for codeml (0.85–0.90) than for Coevol (0.54). We attribute these differences to the autocorrelative or smoothing effect of the Brownian motion model implemented by Coevol [41, 50, 82].

Is ODAM inactivation neutral or adaptive?

The degeneration of morphological structures such as limbs, teeth, and eyes is a complex process that may result from relaxed selection (neutral evolution), adaptive evolution (direct natural selection to conserve energy and/or eliminate the disadvantageous effects of a morphological structure), and/or pleiotropy (indirect selection on another trait) [83]. Adaptive and non-adaptive causes have also been proposed for the inactivation of different genes [28]. In the case of individual lineages of toothless and enamelless placental mammals, it remains unclear if ODAM inactivation resulted from relaxed selection and/or adaptive evolution. For example, was the evolution of edentulism in pangolins driven by the accumulation of random mutations after relaxed selection on tooth development or did natural selection directly favor tooth reduction? If the latter, was this driven by developmental signaling changes that arrested tooth development [84] or were genes that encode structural proteins of the enamel matrix the initial targets of adaptive evolution? This topic is ripe for future studies, but boutique organisms such as pangolins, anteaters, and baleen whales are much less tractable than model organisms such as mice for evo-devo studies of tissue-specific gene expression.

Conclusions

Nature’s laboratory has provided us with multiple independent lineages of placental mammals that are either toothless or lack enamel caps on their teeth. Molecular evolutionary analysis of candidate tooth-specific genes in these independent lineages has emerged as a powerful approach to test hypotheses that various tooth-related genes are pleiotropic versus tooth-specific with respect to their essential functions that are maintained by natural selection. Previous studies have demonstrated that nine different tooth-related genes are inactivated in one or more clades of mammals that are either edentulous or lack enamel caps on their teeth. Here, we show that ODAM is also inactivated in all lineages of toothless and enamelless placental mammals that were investigated. These results support the hypothesis that the only essential functions of ODAM that are maintained by natural selection are related to tooth development and the maintenance of a healthy junctional epithelium where the gingivae are in contact with the tooth enamel. Specifically, ODAM has been inactivated in all representatives of three toothless clades (Vermilingua, Mysticeti, Pholidota) and three enamelless clades (Folivora, Cingulata, Orycteropus) that were included in our study. The overlap in naturally occurring gene inactivations and human genetic diseases [85] suggests that ODAM may be linked to some dental or gingival deformities. DN/dS analyses further demonstrate that the ODAM gene has evolved neutrally in clades that post-date the occurrence of an inactivating mutation on an earlier branch of the placental tree. DN/dS values for transitional branches, which record the first inactivating mutation in each lineage, also provide the basis for estimating ODAM inactivation times. In the main, these estimates are similar to estimates for the inactivation of ENAM in these lineages. Estimates of ENAM and ODAM inactivation in Pholidota both suggest that we may discover stem pholidotans that are older than Eomanis (~ 47 Ma) and lack teeth or at least the enamel on their teeth.

Beyond inactivating mutations in toothless and enamelless mammals, we also discovered mutations in several other groups including odontocetes, four proboscideans, and Callithrix jacchus. In the case of odontocetes, the inactivation of ODAM may be related to a fully aquatic lifestyle where the antimicrobial functions of the junctional epithelium are reduced. An alternate explanation is that the very thin enamel of odontocetes, which is often prismless, may affect the ability of ODAM to contribute to the adhesive properties of the gingival junctional epithelium. The apparent deletion (no BLAST results) of exon 6 of ODAM in three elephantids plus Mammut americanum, but not in two other elephantids, is unexpected and requires some combination of ILS, hybridization, and/or independent deletions of this exon in different proboscideans. However, proboscideans that lack this exon show no other evidence of inactivating mutations and the remaining nine exons have evolved under purifying selection. These results suggest that a version of ODAM that is encoded by nine exons remains functional in these proboscideans.