Influence of Tertiary paleoenvironmental changes on the diversification of South American mammals: a relaxed molecular clock study within xenarthrans

Abstract

Background

Comparative genomic data among organisms allow the reconstruction of their phylogenies and evolutionary time scales. Molecular timings have been recently used to suggest that environmental global change have shaped the evolutionary history of diverse terrestrial organisms. Living xenarthrans (armadillos, anteaters and sloths) constitute an ideal model for studying the influence of past environmental changes on species diversification. Indeed, extant xenarthran species are relicts from an evolutionary radiation enhanced by their isolation in South America during the Tertiary era, a period for which major climate variations and tectonic events are relatively well documented.

Results

We applied a Bayesian approach to three nuclear genes in order to relax the molecular clock assumption while accounting for differences in evolutionary dynamics among genes and incorporating paleontological uncertainties. We obtained a molecular time scale for the evolution of extant xenarthrans and other placental mammals. Divergence time estimates provide substantial evidence for contemporaneous diversification events among independent xenarthran lineages. This correlated pattern of diversification might possibly relate to major environmental changes that occurred in South America during the Cenozoic.

Conclusions

The observed synchronicity between planetary and biological events suggests that global change played a crucial role in shaping the evolutionary history of extant xenarthrans. Our findings open ways to test this hypothesis further in other South American mammalian endemics like hystricognath rodents, platyrrhine primates, and didelphid marsupials.

Background

Paleobiological studies aim at assessing the influence of past environmental changes on the evolutionary history of organisms. To be reliable these macroevolutionary studies must consider groups possessing a high quality fossil record coupled with a well documented history of past environmental changes encompassing geological, ecological, climatic and biogeographical events [1]. Inferring the impact of these historical events on the evolution of faunas is particularly difficult in cases where dispersal occurred repeatedly between biogeographic areas. For this reason, islands have for long attracted particular interest by acting as natural laboratories for the study of evolution [2]. Owing to its "splendid isolation" from other continental land masses for a great part of the Tertiary era [3], South America offers a special case of continental island-like evolution with its very peculiar mammalian fauna composed of an initial nucleus of autochthonous taxa subsequently enriched by few successive waves of successful immigrants.

Numerous studies have attempted to relate the evolutionary history of mammalian fossil faunas in South America with well documented environmental changes [4–6]. However, these studies primarily rely on the quality of the paleontological record which is by essence incomplete. Since the seminal paper of Patterson and Pascual [7] several advances have been made on the understanding of the fossil record from this Southern hemisphere continent [8]. This includes the discovery of additional fossil sites from the previously poorly sampled Tropical zone and new faunal horizons from key ages filling gaps in the stratigraphic sequence. These advances coupled with progress from multidisciplinary studies encompassing tectonic [9], isotopic [10] and radiochemical dating [11] evidences have shed new light on the biotic and environmental history of South America.

Armadillos, anteaters, and sloths belong to Xenarthra and represent the only placental group with living representatives from the initial South American mammalian stock [7]. Owing to their longstanding evolutionary history on this island continent, xenarthrans constitute a promising model for studying the influence of past environmental changes on living organisms. Indeed, the earliest records of Xenarthra in the Paleocene, about 58 million years ago (Mya) [12, 13] already show the morphological features proper to the group (homodonty with simple hypsodont teeth lacking enamel) and possess armadillo-like carapace scutes, suggesting that the group likely originated earlier during the Paleocene or even possibly during the Late Cretaceous [14]. Xenarthra subsequently underwent an impressive Tertiary radiation which led to a great diversity of fossil forms [7]. This diversification was promoted by the South America's isolation from other continental masses at that time. Indeed, South America separated completely from Africa in the Late Cretaceous (80–65 Mya) although it maintained some connections with Antarctica until ca. 36 Mya [15]. Its late Tertiary reconnection with North America via the Isthmus of Panama (3.0–2.5 Mya) led to the Great American Biotic Interchange [16]. During this period of isolation, xenarthrans were one of the dominant groups of placental mammals in the successive South American faunas. They occupied a wide range of ecological niches with fossorial, terrestrial and semiarboreal taxa, gigantic forms such as armoured glyptodonts and ground sloths [17], and sloths even extending their adaptations to the aquatic environment [18]. The late arrival of North American immigrants seems to have only slightly affected xenarthrans and giant terrestrial forms noteworthy became successful southern invaders of Central and North America [16]. Nevertheless, a great part of the xenarthran diversity disappeared with the dramatic extinctions that took place at the end of the Pleistocene [7] and essentially affected the largest forms [19]. This massive and global extinction event left only 30 living xenarthran species whose paleontological origins and phylogenetic relationships are very difficult to decipher from the available fossil record [7], although progresses have recently been made for understanding the origin of some groups of living armadillos [20].

Molecular phylogenetics is now providing a new tool to complement paleontological and geological studies by allowing to reconstruct both the phylogenetic relationships of living species and the timing of their diversification [21]. Indeed, it is likely that global environmental changes left their footprints in the DNA of living organisms from which the evolutionary history can be reliably inferred [22]. Furthermore, there is at present accumulating molecular evidence that past environmental changes have shaped the evolution of diverse terrestrial organisms such as frogs [23], squirrels [24], eulipotyphlan mammals [25] and elephant-shrews [26]. A robust molecular picture of living xenarthran phylogenetic relationships have been recently established [27, 28] and estimates of their divergence dates have already been used to test biogeographic hypotheses [29]. However, the recent development of improved molecular dating methods offers new possibilities in estimating divergence times and allows to go one step forward [21].

Recent methodological advances in molecular dating have aimed at relaxing the assumption of the global molecular clock [30] by modelling the variation of the evolutionary rate along the phylogenetic tree [31–35]. Among these newly developed methods, the Bayesian relaxed molecular clock approach of Thorne et al. [32] appears particularly promising. Indeed, it allows the calculation of divergence time estimates in the presence of rate variation among lineages coupled with the possibility of incorporating multiple paleontological constraints used as priors [36]. This method has already been successfully used to address long standing questions such as the age of the ancestor of the HIV-1 pandemic strain [37] and the diversification of placental mammals in relation to the Cretaceous/Tertiary (K/T) boundary [38–40]. By providing new molecular time scales, the application of the Bayesian relaxed clock has also revealed unsuspected biogeographic patterns in frogs [23] and unravelled the origins of the Malagasy mammalian fauna [41]. These method are nonetheless still in their infancy and robustness of their assumptions to potential sources of errors has to be tested further [34, 42]. However, a significant advance has recently been made with the implementation of the Bayesian method to allow the analysis of genomic data sets by accounting for differences in evolutionary dynamics among genes [43]. Here, we showed that divergence dates within xenarthrans can be accurately inferred by this method and used to unveil the influence of Tertiary environmental changes on the diversification of these South American endemics.

Results and discussion

Evolutionary rate variations and molecular dating methods

The molecular clock hypothesis is rejected for the combined DNA data set by a likelihood ratio test under the F84 + Γ₈ model [44]: δ = 2 x (lnL _{UNCONSTRAINED} - lnL _CLOCK) = 2 x [-87,895.39 - (-89,379.09)] = 2967.40; d.f. = 48; P < 0.0001. This reveals that extensive rate variation occurs among lineages, with fast evolving taxa like Tonatia (Laurasiatheria), murid rodents (Euarchontoglires) and Procavia (Afrotheria), and slow evolving ones like xenarthrans (Xenarthra), perissodactyls (Laurasiatheria), and Dugong and Elephas (Afrotheria) distributed all over the tree. The existence of such an extensive rate variation among lineages precludes the application of the linearized tree method [45] whose estimates can be highly biased when the molecular clock does not hold [46]. Thus, the use of a relaxed molecular clock approach designed to accommodate rate variation is preferable for estimating divergence ages with this data set.

The molecular clock assumption can be relaxed in some ways by defining sub-groups in the tree that evolve locally clocklike and estimating divergence dates within a likelihood framework [47, 48]. Despite a definitive improvement in estimating divergence dates in the case of rate heterogeneity, these methods present some practical limitations pertaining to the arbitrary set up of local clocks and the use of a single calibration point [46, 49]. The management of simultaneous calibration points recently introduced in the local clock method appears nevertheless promising [50]. Another class of methods has attained at modelling the evolutionary rate variation along the tree. The semi-parametric rate smoothing method implements a stochastic model of rate change based on the hypothesis of rate autocorrelation between adjacent branches [31, 35] whereas others attempt to model rate variation across the tree in a Bayesian framework [32–34, 36, 43]. Here, we used a Bayesian relaxed molecular clock method first introduced by Thorne et al. [32] that explicitly accounts for (i) paleontological uncertainty on calibration points by estimating the posterior distribution of rates and divergence times from prior distributions on rates and paleontological constraints [36], and (ii) contrasted patterns of evolutionary rate among genomic loci [43] as this has been shown to be the case for our data set [27].

Basal placental mammal divergences and the K/T boundary

Our dataset comprising 47 placental mammals allowed us to estimate divergence dates for the major placental clades in relation to the Cretaceous/Tertiary boundary (K/T). Indeed, there has been considerable debate [51, 52] on the impact of this massive extinction event on early placental mammal lineages [53]. A recent molecular study by Springer et al. [40] using the same Bayesian dating method and calibration points – but for a longer data set – found support for the Long Fuse model of Archibald and Deutschman [54] by placing almost all interordinal divergences in the Late Cretaceous. The results we obtained in estimating the age of the Most Recent Common Ancestors (MRCAs) of seven basal placental clades using different data partitions of our nuclear data set are recapitulated in Table 1. It is noteworthy than the partitions dominated by 1^st and 2^nd codon positions (ABV 1+2 and A2BV2) generally yielded slightly more recent dates than the complete data set, a behaviour also noticed by Springer et al. [40] on nuclear exons. Nevertheless, irrespective of the data partition, all estimates obtained for the MRCAs of basal placental clades fell within the Late Cretaceous (Table 1). As there is only slight differences between the results obtained from the three different partition schemes, we will here focus on the results obtained from the longest dataset.

Table 1 Mean Bayesian posterior estimates for basal placental and xenarthran divergence dates, and calibration points using different gene partitions. Molecular datings were calculated with MULTIDIVTIME 43 using the following gene partitions of the nuclear genes (ADRA2B, BRCA1, and VWF): ABV (1+2) = 1^st and 2^nd codon positions of the three genes (3421 sites); A2BV2 = 1 ^st and 2^nd codon positions of VWF and ADRA2B plus all codon positions of BRCA1 (4350 sites); and ABV (1+2+3) = all codon positions of the three genes (5130 sites). Standard deviations (SD) and 95% credibility intervals (CredI) are indicated. Nodes are named and numbered as in Figures.

The chronogram obtained on the complete placental data set is presented in Figure 1. Based on this data set, the mean posterior age for the placental crown group, corresponding to the early emergence of Afrotheria, is estimated to be in the Early Cretaceous at 107 ± 6 Mya (95% Credibility Interval: 97–119 Mya). This early split was quickly followed by the divergence between Xenarthra and the remaining placentals at the end of the Early Cretaceous 103 ± 5 Mya (94–113). All four major placental clades have their last common ancestors in the Late Cretaceous (Figure 1) at respectively: 82 ± 6 Mya for Afrotheria (71–95), 65 ± 5 Mya for Xenarthra (55–75), 89 ± 4 Mya for Euarchontoglires (81–100), and 83 ± 3 Mya for Laurasiatheria (78–90). The divergence between Euarchontoglires and Laurasiatheria is estimated to have occurred at 94 ± 4 Mya (87–103). On the whole, our mean posterior estimates of basal placental divergence dates are fully compatible with those obtained by Springer et al. [40] on a similar set of characters (i.e. nuclear exons) and using the same priors on calibration points. The main placental divergences appear to have occurred in the Cretaceous (Figure 1) as previously proposed by other studies [55, 56] based on the results from the quartet dating approach [47] or by using the Bayesian relaxed molecular clock approach on complete mitochondrial genomes [38, 39] and concatenated nuclear DNA sequences [25, 40]. The good agreement between these studies in obtaining Cretaceous dates for basal placental divergences seems to validate the occurrence of a massive survival of modern placental lineages across the K/T boundary [40]. However, more paleontological data on the degree of morphological and ecological specialization of early placental lineages are required to distinguish further between alternative models of diversification [53]. In that respect, the recent discovery of the very well conserved eutherian fossil Eomaia in the cretaceous of China (125 Mya) is encouraging [57].

Figure 1

Molecular dating of basal placental divergences and the Cretaceous/Tertiary boundary. Chronogram with branch lengths proportional to time units obtained from the partitioned analysis of the complete ADRA2B + BRCA1 + VWF (codon positions 1+2+3) placental nucleotide dataset (5130 sites). The geological timescale for the Cretaceous and Tertiary eras is superimposed (Paleo. = Paleocene; Oligo. = Oligocene; P. = Pliocene). The thick horizontal dashed line marks the Cretaceous/Tertiary boundary (K/T) at 65 Mya. The grey shaded area representing the Late Cretaceous era illustrates the fact that the major placental divergences occurred within this time window. The names of basal placental clades (yellow diamonds) and calibration nodes (red stars) referring to Table 1 are indicated (PAE: PAEnungulata, MUR: MURidae, MEG: MEGachiroptera, PER: PERissodactyla, ART: CetARTiodactyla, and CET: CETacea). The xenarthran subtree is highlighted in bold and numbering of nodes refers to Table 1.

Another insight from our results is the fact that the paleontological constraint put on the paenungulate MRCA (Figure 1) appears to represent a large underestimate of the actual divergence time for this group (Table 1). In fact, the mean posterior estimate for this node a priori constrained at a minimum of 54 Mya was 69 ± 6 Mya (59–81). This suggests that the earliest proboscidean fossil [58] is far from representing the earliest Paenungulate and that older fossils belonging to this group remain to be found. This situation might be extended to the African clade (Afrotheria) as a whole since its fossil record remains poorly known despite the recent discovery of a spectacular Oligocene fauna [59]. The assessment of the relative reliability of independent paleontological constraints represents a point that certainly needs more scrutiny in molecular clock analyses since it might have a strong impact on the results [46]. The use of the more precise and reliable calibration constraints is indeed advocated to be of outmost importance in the application of Bayesian molecular clock methods [43]. By allowing to reciprocally estimate the posterior age of different nodes for which paleontological constraint are available, the Bayesian molecular clock method might help to identify unreliable calibration points and groups for which older fossils are likely to be found.

A molecular timescale for xenarthran evolution

The age of the xenarthran crown group estimated by the Bayesian method in the present study (65 ± 5 Mya; 55–75) is consistent with previous results obtained from a local molecular clock approach [48] on the VWF alone suggesting a 59 to 76 million years (Myr) interval [29]. Xenarthra thus appears to be the major placental clade with the youngest MRCA, but also with the longest ancestral lineage (Figure 1). The first appearance of xenarthrans in the fossil record is materialized by the occurrence of the earliest armadillo scutes in the late Paleocene (ca. 58 Mya) of Brazil [13, 15]. Therefore, our results are compatible with the age assumed for the xenarthran diversification, but suggest the existence of a xenarthran "ghost" lineage that left no fossil traces for almost 50 Myr since their purported origin estimated here around 105 Mya. This observation concurs with the fact that the origin of xenarthrans still constitutes a paleontological and biogeographic enigma [16, 60]. Only the discovery of new fossils – presumably of Cretaceous age according to our estimates – might help to solve this mystery.

The results obtained for the timing of the xenarthran diversification by the relaxed molecular clock approach are recapitulated in Table 1. Following the early split between Cingulata (armadillos) and Pilosa around 65 Mya at the end of the Cretaceous, anteaters and sloths separated at the transition between Paleocene and Eocene some 55 Mya (46–65). Within anteaters the lineage leading to the pygmy anteater (Cyclopes) emerged in the Middle Eocene around 40 Mya (32–49), whereas the Tamandua and Myrmecophaga lineages diverged 30 Myr later in the Late Miocene at 10 Mya (7–14). The oldest undoubted fossil of anteaters comes from the Colhuehuapian (ca. 20 Mya) South American Land Mammal Age (SALMA) and they are well known since the Santacrucian SALMA (ca. 16 Mya). Our estimates, matching those obtained in Delsuc et al. [29], suggest a long evolution of anteaters prior to the Colhuehuapian SALMA, and confirm the incompleteness of the early fossil record of anteaters. The molecular estimates also emphasize the relative antiquity of the pygmy anteater lineage whose only living representative (Cyclopes didactylus) is morphologically very divergent from the other two genera [61]. The considerable evolutionary divergence of the pygmy anteater might be more efficiently reflected in the taxonomy by placing it in its own family as advocated by Barros et al. [62].

The Bayesian method estimates the separation between the two modern sloth lineages at 21 Mya (15–28). Unknown as fossils [7], the extant genera of two-toed (Choloepus) and three-toed sloths (Bradypus) were once placed into two distinct families (respectively Megalonychidae and Bradypodidae) on the basis of their numerous morphological differences and a presumably diphyletic origin [63]. Recent cladistic analyses based on craniodental evidence, including the living forms and the three traditionally recognized families of fossil sloths (Megatheriidae, Megalonychidae and Mylodontidae), placed Bradypus as the sister-taxon of all remaining sloths and Choloepus with the extinct Megalonychidae [64]. The oldest known sloths come from the Eocene of Patagonia and Antarctica but they cannot be precisely assigned to any of the recognized lineages [65]. Our fairly ancient estimation confirms the considerable divergence between the two modern sloth genera bringing support for a taxonomic distinction at a high rank.

The almost complete generic sampling among armadillos (all living genera but the fairy armadillos Chlamyphorus) allows for the first time an adequate estimation of the divergence dates within this poorly known group of placental mammals. The early emergence of Dasypodinae – including the living long-nosed armadillos (Dasypodini) and the fossil Stegotheriini [13, 66] – within armadillos is estimated to have occurred during the Middle Eocene around 40 Mya (31–49). If the plesiomorphic Astegotheriini are removed from the Dasypodinae (see [13]), this estimation is consistent with the paleontological evidence that suggests an early divergence of this group from other lineages containing living representatives [66]. The two divergent living species of long-nosed armadillos studied (Dasypus novemcinctus and D. kappleri) would have separated at 7 Mya (5–11) in the Late Miocene. This is an interesting result that obliges paleontologists to review the fossil evidence and could also be significant in the taxonomy of the living species. The earliest and most primitive Dasypodini (Anadasypus hondanus) comes from the middle Miocene Laventan SALMA (ca. 13 Myr) in Colombia [20]. That means that the origin and radiation of the Dasypus genus must have occurred very soon after the emergence of the group and paleontologists should expect to find Dasypus remains at any stratigraphic level younger than the Late Miocene. However, the hard evidence of Pre-Pleistocene record is very weak. The oldest undoubted fossil record of the genus comes from the Late Pleistocene, around 0.8 Mya [67]. This fossil evidence would suggest that the separation between D. novemcincus and D. kappleri should be much younger than proposed here, during the Pleistocene or even the Holocene. However, this can be very easily a matter of the paleontological record. Although the early existence of Dasypodinae is recorded in Patagonia, the origin and evolution of the Dasypodini seems to have happened mostly in Northern South America [20], where the paleontological exposures of appropriate age are not as spectacular as in Southern South America.

The split between the two remaining armadillo sub-families Tolypeutinae (Tolypeutes, Cabassous, and Priodontes) and Euphractinae (Euphractus, Chaetophractus, and Zaedyus) took place about 33 Mya (25–42). The respective diversifications of these two subfamilies appear to have happened in relatively short time periods between 20–22 Mya for Tolypeutinae and 6–7 Mya for Euphractinae (Table 1; Figures 2,3). The evolution of the Tolypeutinae seems to parallel that of the Dasypodini [20]. This is particularly evident for the Tolypeutini, including the earliest known Pedrolypeutes praecursor from the Middle Miocene Laventan SALMA in Colombia, and the genus Tolypeutes, known since the Pleistocene of the pampean region. The fossil record of the giant armadillo (Priodontes) and naked-tailed armadillos (Cabassous) is virtually unknown, but the present distribution clearly suggests a Northern South American history as for the Tolypeutini and Dasypodini. On the other hand, the evolution of the euphractines is more related to Southern South America. The early history of Euphractini is recorded in Patagonia since the Deseadan SALMA (30–25 Mya) with Prozaedyus a fossil form that clearly resembles the living members of the group. However, Carlini and Scillato-Yané [68] proposed that the Late Miocene-Pliocene euphractines represent an independent radiation in temperate to warm environments of central and northern Argentina from that of the living genera.

Figure 2

Synchronicity of diversification events among independent xenarthran lineages. The histograms represent the distributions of 10,000 posterior divergence time estimates sampled during the MCMC analysis of the complete ADRA2B + BRCA1 + VWF (codon positions 1+2+3) placental dataset (5130 sites). A. First split within anteaters (node 9) versus first split within armadillos (node 10). The smaller histogram (top right) shows the posterior distribution of the difference between the ages estimated for these nodes (node 10 – node 9). B. Split between three-toed and two-toed sloths (node 6) versus diversification of tolypeutine armadillos (nodes 5 and 7). C. Split between the two divergent species of long-nosed armadillos (node 3) versus diversification of euphractine armadillos (nodes 1 and 2). Numbering of nodes refers to Table 1 and Figures 1 and 3.

Figure 3

Molecular timing of the xenarthran radiation and Tertiary major environmental changes in South America. Phylogenetic relationships and divergence times are represented by a chronogram, whose branch lengths are proportional to time units, obtained from the partitioned analysis of the complete ADRA2B + BRCA1 + VWF (codon positions 1+2+3) placental dataset (5130 sites). Confidence limits on divergence time estimates are represented by rectangles at nodes corresponding to ± one standard deviation (see Table 1). Tertiary Epoch boundaries follow the 1999 geologic timescale of the Geological Society of America [77]. The Cretaceous/Tertiary transition (K/T) is represented by a vertical dashed bar at 65 Mya. Major Andean tectonic crises 5 are represented by grey shaded areas. Eustatic curves of sea level are from [71]. The ocean temperature curve (in red) is based on the high-resolution deep-sea oxygen isotope record (δO¹⁸) [72]. Node numbering refers to Table 1.

These results confirm the fact that living Dasypodidae contains quite divergent taxa grouped into subfamilies of fairly ancient origins [29]. This suggests that the still unclear paleontological roots of these groups have to be searched deeper than previously thought [7]. The relatively younger dates obtained within armadillos in Delsuc et al. [29] might be explained by an effect of reduced taxon sampling since only three genera (Dasypus, Cabassous and Chaetophractus) were previously considered. The ancient dates obtained for the origin of each of the three armadillo subfamilies are coherent with their distinctive morphologies [69] and marked differences observed in the structure of their spermatozoa [70].

Xenarthran diversification and paleoenvironnemental changes

In the light of these results, it is striking to note that some diversification events within Xenarthra seem to be synchronous. This synchronicity among independent xenarthran lineages is illustrated by the correspondence between the distributions of posterior divergence time estimates sampled during Markov Chain Monte Carlo (MCMC) runs presented in Figure 2. Indeed, the separation of the Cyclopes lineage from other anteaters is perfectly correlated with the separation of the Dasypus lineage from other armadillos around 40 Mya in the Middle Eocene. The quasi perfect correspondence between these two speciation events is well illustrated by looking at the posterior distribution of the differences between the ages estimated for these nodes which is centred around 0 (Figure 2A). Similarly, the diversification of the armadillo subfamily Tolypeutinae (Priodontes, Cabassous and Tolypeutes) relates to the separation between the two modern sloths genera around 21–22 Mya in the Early Miocene (Figure 2B). And finally, the recent diversification of the subfamily Euphractinae (Chaetophractus, Euphractus and Zaedyus) also matches well with the separation between the two living species of the genus Dasypus here included at 6–7 Mya in the late Miocene (Figure 2C). Such a correlated history between independent lineages is unlikely to have occurred by chance alone and suggests a role for external factors in shaping the evolutionary history of living xenarthrans.

The influence of environmental changes that occurred during the Cenozoic of South America on the evolution of its endemic fauna has long been studied from the paleobiogeographical point of view [7]. Indeed, this era was characterized by drastic climatic variations associated with changes in sea levels [71], major biotic modifications [4] and tectonic phases of Andean uplift [5]. Moreover, recent spectacular advances have greatly improved the understanding of the fossil record of South America [8]. These advances from multidisciplinary studies on tectonics [9], stable isotopes [10] and biostratigraphy and geochronology [11] have clarified the biotic and environmental history of South America during the Cenozoic. The making up of this integrated framework has rendered possible to retrace and identify general evolutionary patterns for terrestrial communities such as mammalian herbivores [6]. Until recently, there has been only limited efforts to incorporate the results of molecular studies into this multidisciplinary canvas [22]. However, some recent surveys have provided substantial evidence that past environmental changes may have significantly influenced the evolutionary history of mammalian groups at a large scale [24–26].

As the only South American endemic group of placental mammals whose evolutionary history encompasses all the Cenozoic, Xenarthra represents a good candidate to test the hypothesis that global environmental changes have influenced the diversification of the endemic fauna. Interpreting our results in the context of relatively well documented environmental variables like Andean tectonic crises, changes in ocean level, and temperature variations sheds light on the peculiar diversification pattern revealed by our molecular dating analyses. Indeed, each of the three correlated diversification events identified among independent xenarthran lineages appears to follow periods of major environmental changes (Figure 3). More specifically, these events seem to be slightly posterior to the major phases of Andean uplift defined as "tectonic crises" by Marshall and Sempere [5]. First, the codiversification of Myrmecophagidae and Dasypodidae in the Middle Eocene occurred just after the large Incaic uplift episode (Figure 3). The estimated date of ca. 40 Mya for these synchronous diversification events matches up with a well dated pulse of this deformation episode in the Andes of Peru estimated at ca. 43 Mya [5]. Second, the simultaneous diversification of modern sloth lineages and the proposed radiation of Tolypeutinae in the middle of the Early Miocene correlates with the end of the first major Bolivian crisis (Figure 3). This diastrophic event that occurred during the Late Oligocene-Early Miocene interval was an intense deformational and magmatic episode widespread along the Andes [9]. The intensity of this crisis is seen as a turning point in Andean tectonics because the Andes became at that time the principal relief of the West coast of South America creating a rain shadow that significantly influenced South American climates [5]. This period is also notably marked by a thermal optimum followed by a brief but deep glacial maximum [72] and sharp marine regressions [71]. This leads to a major shift in South American mammalian fossil communities with the beginning of the Miocene radiation of ground sloths for example [7]. Third, both the radiation of living euphractine armadillos and the separation between long-nosed armadillos correspond with the latest phases of the second major Bolivian tectonic crisis marked by a significant increase in Andes height (Figure 3). This period marks a stage of general cooling with the formation of the Antarctic ice sheet [72] resulting in a global drying of the Patagonian region and the formation of the Argentinean pampas in the Late Miocene [7]. This major shift in vegetation with the beginning of the predominance of savannas and grasslands was recorded in isotopic studies in the form of an elevation of the C₄/C₃ ratio [10]. The fact that euphractine armadillos represent a zoogeographically well defined group that occupies the pampas and savannas of Southern South America [69] adds credit to the idea that this major environmental change created ecological opportunities that triggered the diversification of this subfamily. The evolutionary history of extant xenarthrans also appears to have been influenced by the major environmental changes that marked the transition between the Eocene and Oligocene epochs with the split at this time between Tolypeutine and Euphractine armadillos (Figure 3). This period corresponds to a glacial optimum with the first formation of a concrete Antarctic ice sheet and the creation of the circum-Antarctic oceanic current that changed the atmospheric circulation and induced a drastic general cooling [72]. This marked the passage from warm, humid, tropical-temperate forest environments to much more arid and dryer habitats dominated by savannas and grasslands [4] prompting the diversification of mammalian herbivore communities [6].

Overall, our results are consistent with the hypothesis of a major influence of global environmental change on the evolution of xenarthrans. This process appears to have been apparently mainly driven by the consequences of Andean tectonics as suggested by Marshall and Sempere [5] even if the causal link between climate change and the rise of the Andes might be more complicated than previously though [73]. These authors argued as a working hypothesis that geological and tectonic processes might have played a key role in restructuring the mammalian communities by drastically remodelling the South American environments and climates. Some groups have consequently been driven to extinction and new ecological opportunities have been provided for others to diversify. In particular, the finding that the two radiation events identified within armadillos [28] appear to follow two major tectonic crises is consistent with Marshall and Sempere's hypothesis and suggests that other living members of South American endemic groups might also have recorded these events in their genomes.

Conclusions

We have used a Bayesian relaxed molecular clock approach explicitly taking into account paleontological uncertainty and contrasted evolutionary dynamics between genomic loci to obtain divergence time estimates for living xenarthrans. We proposed a time scale for the diversification of this major placental clade with a sparse fossil record, based on the analysis of three nuclear genes for 50 mammals. This molecular timescale put in relation with well documented environmental changes that occurred during the Cenozoic of South America revealed the crucial influence that global change – possibly induced by Andean tectonic processes – might have played on the evolutionary history of these peculiar mammals. It is finally suggested that these major paleoenvironmental changes may also have left their footprints in the genomes of other South American endemic mammalian taxa. This hypothesis has yet to be tested in hystricognath rodents, platyrrhine monkeys or didelphid marsupials for example.

Methods

Data sets and topology

The original data set we use here is the one obtained by Delsuc et al. [27] consisting of 50 mammalian taxa for the three nuclear genes α2B adrenergic receptor (ADRA2B), breast cancer susceptibility exon 11 (BRCA1), and von Willebrand factor exon 28 (VWF). The complete alignment of these genes represents a total of 5130 nucleotide sites [27]. Different partitions from this data set were used for calculating divergence times: ABV (1+2) = concatenation of first and second codon positions of the three genes (3421 sites); A2BV2 = concatenation of first and second codon positions of ADRA2B and VWF plus all codon positions of BRCA1 (4350 sites); and ABV (1+2+3) = concatenation of all codon positions of the three genes (5130 sites). As the relaxed molecular clock approach relies on a topology to infer divergence times, we used the ML topology previously identified on the complete data set [27] which conforms to the current views on placental mammals relationships based on the latest large scale analyses [74].

Calibration points

In the Bayesian relaxed molecular clock approach [32], it is important to use prior constraints on independent calibration points dispersed across the tree in order to reduce potential regional effects [36]. To render things comparable among studies, we used the calibrations of Springer et al. [40] that are compatible with our taxon sampling. Thus, the six calibration points defining the following eight prior constraints were used (see Figure 1): (1) a minimum of 54 Mya for Paenungulata (PAE), (2) a minimum of 12 Mya for the divergence between the Muridae Mus and Rattus (MUR); (3) a minimum of 43 Mya and a maximum of 60 Mya for Megachiroptera (MEG); (4) a minimum of 54 Mya and maximum of 58 Mya for Perissodactyla (PER); (5) a conservative maximum of 65 Mya for Cetartiodactyla (ART); and (6) a minimum of 52 Mya for the first appearance of Cetacea (CET). We do not used prior constraint on the origin of Xenarthra because we wanted to make independent inferences about their divergence dates. We verified however that constraining the base of Xenarthra with a minimum of 60 Mya [40] does not significantly affect our molecular estimates (data not shown). For the same reasons of independence, we did not use intra Xenarthra calibrations, but also because of a lack thereof. Indeed, if the fossil record of xenarthrans is especially rich for extinct giant sloths for example, the paleontological origins of extant species are almost unknown [7]. In all subsequent Bayesian analyses these eight prior constraints on calibration points were simultaneously used to derive posterior estimates of divergence ages [36].

Divergence time estimates

The Bayesian relaxed molecular clock approach was applied using the program package MULTIDIVTIME [75]. The module ESTBRANCHES was first used to estimate branch lengths of the constrained topology and the corresponding variance-covariance matrices for each data partition using the three marsupials Macropus, Didelphis and Vombatus as outgroups for the 47 placental taxa. The F84+Γ₈ model was used with maximum likelihood parameters previously estimated by PAML version 3.13 [76]. Then, the module MULTIDIVTIME used the variance-covariance matrices produced by ESTBRANCHES to run a MCMC for estimating mean posterior divergence times on nodes with associated standard-deviation (SD) and 95% credibility interval (CredI). In all calculations, the MCMC was sampled 10,000 times every 100 cycles and the burn-in stage was set to 100,000 cycles in all subsequent analyses. The following priors were used for the mammalian data set: 110 Mya (SD = 55 Mya) for the expected number of time units between tip and root if there has been no constraint on node times, and 200 Mya for the highest possible number of time units between tip and root. Other priors for gamma distribution of the rate at root node and the Brownian motion constant describing the rate variation (i.e., the degree of rate autocorrelation along the descending branches of the tree), were derived from the median branch length for each data set as advised by Thorne et al. [32].