Evolutionary Relationships among Extinct and Extant Sloths: The Evidence of Mitogenomes and Retroviruses.

Macroevolutionary trends exhibited by retroviruses are complex and not entirely understood. The sloth endogenized foamy-like retrovirus (SloEFV), which demonstrates incongruence in virus-host evolution among extant sloths (Order Folivora), has not been investigated heretofore in any extinct sloth lineages and its premodern history within folivorans is therefore unknown. Determining retroviral coevolutionary trends requires a robust phylogeny of the viral host, but the highly reduced modern sloth fauna (6 species in 2 genera) does not adequately represent what was once a highly diversified clade (∼100 genera) of placental mammals. At present, the amount of molecular data available for extinct sloth taxa is limited, and analytical results based on these data tend to conflict with phylogenetic inferences made on the basis of morphological studies. To augment the molecular data set, we applied hybridization capture and next-generation Illumina sequencing to two extinct and three extant sloth species to retrieve full mitochondrial genomes (mitogenomes) from the hosts and the polymerase gene of SloEFV. The results produced a fully resolved and well-supported phylogeny that supports dividing crown families into two major clades: 1) The three-toed sloth, Bradypus, and Nothrotheriidae and 2) Megalonychidae, including the two-toed sloth, Choloepus, and Mylodontidae. Our calibrated time tree indicates that the Miocene epoch (23.5 Ma), particularly its earlier part, was an important interval for folivoran diversification. Both extant and extinct sloths demonstrate multiple complex invasions of SloEFV into the ancestral sloth germline followed by subsequent introgressions across different sloth lineages. Thus, sloth mitogenome and SloEFV evolution occurred separately and in parallel among sloths.

Introduction

Superorder Xenarthra comprises one of the most remarkable of the major extant groups of placental mammals. Its constituent orders include Folivora (sloths) and Vermilingua (anteaters), which together form the sister group to Cingulata (armadillos and relatives). Although still an important part of the South American mammal fauna, much about their phylogenetic history remains profoundly obscure (Murphy et al. 2001; Meredith et al. 2011; Dos Reis et al. 2012; O’Leary et al. 2013). Endogenous retroviruses (ERVs) represent ancient viral infections that have become integrated into a host germline and are thereafter transmitted vertically as simple Mendelian traits. Widespread ERVs within a taxon may be considered to be fixed. ERVs evolve neutrally (Kijima and Innan 2009; Feschotte and Gilbert 2012), at rates much slower than those of their exogenous counterparts. Optimally, rates exhibited by retroviral loci and those of the host germline will be in essential synchrony, thus permitting their potential use as codivergence markers for phylogenetic studies. A group of ERVs known as foamy viruses are widely distributed among eutherian mammals and appear to have undergone coevolution and cospeciation in synchrony with their hosts for more than 100 Myr (Katzourakis et al. 2014). Foamy viruses are nonpathogenic, complex retroviruses and form a unique group among the Retroviridae, Spumaretrovirinae (Linial 1999).

A recent study of ERVs in extant sloths concluded that sloth endogenized foamy virus (SloEFV) invaded the genome of a sloth ancestor in the late Middle Eocene some 39 Ma (Katzourakis et al. 2009), before the estimated divergence of two- and three-toed sloths (∼21 Ma) but subsequent to the anteater-sloth divergence (∼55 Ma). This inference is of interest because sloths comprise a formerly diverse group of placental mammals, with nearly 100 species described in the fossil record (McKenna and Bell 1997; fig. 1). Although clearly an ancient group, they were mostly limited to South America until mid- to late Cenozoic times, when certain lineages reached various islands in the West Indies and North America (Steadman et al. 2005; MacPhee et al. 2007; Gaudin and McDonald 2008; McDonald and De Iuliis 2008). Paleontologically, the oldest reasonably well preserved fossils assigned to Folivora are Oligocene (Carlini and Scillato-Yané 2004; Pujos and De Iuliis 2007), although it must be emphasized that better support is needed to test whether these earliest sloths actually fit cladistically among the crown group (see Discussion). Declining post-Miocene diversity eventually culminated during the late Quaternary in the loss of all sloth taxa, except for the surviving tree sloths: The three-toed Bradypus (Bradypus tridactylus, Bradypus variegatus, Bradypus torquatus, and Bradypus pygmaeus) and the two-toed Choloepus (Choloepus hoffmanni and Choloepus didactylus) sloth species.

F

Constraint tree used for divergence time estimation in BEAST under the fossilized birth-death process. Extinct taxa are depicted at the midpoints of their stratigraphic ranges (supplementary table S2, Supplementary Material online) though analyses integrate over uncertainty regarding true ages. Partitioning of the tree into named taxonomic groups at the family and higher levels is provided for general reference.

Given how restricted the extant complement of sloths is to make any headway in understanding the evolution of ERVs within folivorans, it would be highly desirable to investigate their incidence in extinct forms. Initial attempts to obtain ancient DNA (aDNA) sequences from different species of Late Pleistocene terrestrial (“ground”) sloths for phylogenetic purposes were broadly successful (Höss et al. 1996; Poinar et al. 1998, 2003; Greenwood et al. 2000; Hofreiter et al. 2003; Clack et al. 2012), suggesting that it may be possible to obtain sequence information on SloEFVs occurring in extinct species. A barrier to better understanding of SloEFV evolution within this group remains because morphological (Gaudin 1995, 2004; White and MacPhee 2001) and previous molecular analyses of sloth phylogeny have failed to reach a consensus on higher level relationships among extinct and extant taxa.

To clarify phylogenetic relationships and estimate divergence dates among and between sloth lineages, as well as to determine SloEFV macroevolutionary patterns within folivorans, we determined complete mitogenomes and a portion of the SloEFV pol gene using a hybridization capture enrichment technique (Maricic et al. 2010) and high-throughput sequencing. This methodology has been successfully used to retrieve full mitochondrial genomes and nuclear loci from samples ranging from modern to tens of thousands of years old (Tsangaras et al. 2014; Sarkissian et al. 2015). Our results corroborate previous molecular estimates of higher level folivoran relationships but suggest a more complex co-volutionary history between SloEFVs and hosts than has been previously suggested.

Materials and Methods

Preliminary Issues

The endogenization of ERVs is a complicated process in which an exogenous retrovirus initially infects a host and then undergoes continuous amplification, reinfection and recolonization under host selection pressure before either being removed from the population by drift or becoming fixed in the host genome as an ERV (Gifford and Tristem 2003). However, there are several potential problems with the use of fixed ERVs for phylogenetic purposes. It is possible, for example, that different insertions of the same ERV found in a host may have originated from multiple independent infections or introgression events as ERV-containing individuals breed with ERV-free individuals. Such a process of introgression is seen in the case of the retrovirus (KoRV), currently undergoing endogenization in koalas (Phascolarctos cinereus) (Ishida et al. 2015). Evidence of host switching and introgression have been found among distantly related vertebrates (Hayward et al. 2013, 2015) as well as among closely related mammals (Switzer et al. 2005; Jern et al. 2006; Katzourakis et al. 2014). In consequence, phylogenetic dissociation between host and ERVs certainly occurs, but detecting it from patterns of infection or introgression in extant taxa can be difficult because of the long time spans involved and because many hosts and exogenous retroviral counterparts within taxonomic groups are extinct. However, when retroviral and host evolution are synchronous, retroviruses can be useful in determining host phylogenetic relationships to the study of cross-species transmission and retroviral introgression patterns require a well-resolved host phylogeny.

Samples

An ∼13,000-year-old Mylodon darwinii bone specimen (previously described in Greenwood et al. 2000) and a Mylodon coprolite sample (age not determined) were obtained from the Natural History Museum, London (sample number BM(NH)M8758) (table 1). Two ∼20,000-year-old coprolite samples attributed to Nothrotheriops shastensis, one from Gypsum (USNM 617523), Nevada, and the other from Rampart Cave (USNM 617524), Arizona, were provided by the National Museum of Natural History, Smithsonian Institution. Tissue (spleen) of C. didactylus was obtained from a sloth from the Tierpark Berlin which was dissected post mortem at the Leibniz Institute for Zoo and Wildlife Research, Berlin. Fresh fecal samples from two- and three-toed sloths (C. didactylus, C. hoffmanni, and B. variegatus, respectively) were obtained in the wild at Playa Bluff, Bocas Del Toro, Panama. Samples were collected and stored in alcohol and frozen at −20 °C until processing.

Table 1

Sloth Sample Information

Species	Sampling Location/Provider	Description	Experimental Use
Choloepus didactylus	Tierpark Berlin, collected 2006	Spleen	SloEFV pol gene and sloth mitochondrial genome
C. didactylus	Tierpark Berlin, collected 2014	Multiple tissues	SloEFV pol gene
C. didactylus	Tierpark Berlin, collected 2010	Coprolite	Hybridization capture bait generation
Choloepus hoffmanni	Smithsonian Tropical Research Institute, collected in Bocas del Toro, Panama	Coprolite	SloEFV pol gene and sloth mitochondrial genome
Bradypus variegatus	Smithsonian Tropical Research Institute, collected in Bocas del Toro, Panama	Coprolite	SloEFV pol gene and sloth mitochondrial genome
Mylodon darwiniia	BM(NH)M8758, Natural History Museum, London	Bone	SloEFV pol gene and sloth mitochondrial genome
M. darwinii	Natural History Museum, London	Coprolite	SloEFV pol gene and sloth mitochondrial genome
Nothrotheriops shastensis	USNM 617523, National Museum of Natural History, Smithsonian Institution, Quaternary/Pleistocene Gypsum Cave, Nevada	Coprolites	SloEFV pol gene and sloth mitochondrial genome
N. shastensis	USNM 617524, National Museum of Natural History, Smithsonian Institution, Quaternary/Pleistocene Rampart Cave, Arizona	Coprolites	SloEFV pol gene and sloth mitochondrial genome

The same sample is used by Greenwood et al. (2000) and is estimated to be 13,000 year old.

DNA Extractions

DNA extractions from extinct sloths were performed in a facility dedicated to historical DNA and aDNA extractions in the Department of Wildlife Diseases of the Leibniz Institute for Zoo and Wildlife Research under a laminar flow hood and ultraviolet (UV) decontamination. The laboratory has never been used for molecular work on modern genetic samples. The room is UV irradiated 4 h every night from ceiling mounted UV lights. Protective clothing was always worn to avoid contamination derived from personnel. From the bone material, powder for extraction was obtained by drilling with a low speed hand-held drill. Coprolite samples were ground to powder with a Freezer/mill 6700 bone grinder.

Approximately 0.25 g of bone or coprolite powder was extracted using a silica-based extraction kit for aDNA (GENECLEAN Ancient DNA Extraction Kit; MP Biomedicals). The protocol followed the vendor’s instructions and has been successfully applied to a variety of ancient sample types (Wyatt et al. 2008; Roca et al. 2009). Mock extractions were performed with each museum specimen as negative control for contamination during extraction. DNA extracts of each museum specimen was further purified using MinElute spin columns (Qiagen, Hilden, Germany) as described in Gilbert et al. (2007) to remove potential inhibitors for the subsequent enzymatic reactions.

For modern samples, total DNA was extracted from ∼25 mg of spleen, of C. didactylus using the Qiagen DNeasy Blood & Tissue Kit (Qiagen) following the manufacturer’s protocol. For each of the modern coprolite samples, 200 mg was used for extraction using the PSP Spin Stool DNA Kit (STRATEC Biomedical, Birkenfeld, Germany). All extraction work on modern samples was performed under a laminar flow hood in a laboratory separate from that used for the ancient samples.

Ethics Statement

Experiments involving sloth tissues were approved by the Internal Ethics Committee of the Leibniz Institute for Zoo and Wildlife Research, approval number 01-01-2013. The fresh fecal samples for modern sloth from Panama were collected under ANAM permit SE/A-61-10.

Preparation of Baits and Amplicons

The mitochondrial genome hybridization capture bait was generated using DNA from a spleen extract of C. didactylus. Four pairs of primers were designed to cover the sloth mitochondrial genome in four overlapping fragments using Primer3web version 4.0.0 (Sl mt F1: 5′-CAGCCCATGATCACACATAAC-3′; Sl mt R1: 5′-TCCTCAGCCYCCAAYTATAA-3′; Sl mt F2: 5′-CAACAGAAGCAGCCACCAAATA-3′; Sl mt R2: 5′-GTTRAYTATGTGGTA GGCGTGTA-3′; Sl mt F3: 5′-AARACCAAARCATCACTAGCC-3′; Sl mt R3: 5′-CATCCG ATYAGTAGGAARGAT-3′; Sl mt F4: 5′-ATTCAYCCCAGTAGCACTC-3′; Sl mt R4: 5′-GTCGATTATAGGACAGGTTCC-3′). Polymerase chain reactions (PCRs) were performed using Platinum® PCR SuperMix High Fidelity (Life Technologies, Darmstadt, Germany) with protocol including an initial 94 °C for 4 min; 35 cycles at 94 °C for 30 s, 54 °C for 30 s, and 72 °C for 6 min; final extension of 8 min at 72 °C. An aliquot of each PCR product was visualized on a 1% agarose gel stained with GelRed (Biotium, Dossenheim, Germany). PCR products were purified using QIAquick columns (Qiagen) and then Sanger sequenced to verify the products. The verified PCR products were sheared using a microTUBE (Covaris Ltd., Brighton, UK) to an average length of 150 bp. The sheared DNA was quantified using a NanoDrop ND-1000 (Thermo-Scientific) and the equimolar products pooled.

SloEFV pol gene sequences were generated using the DNA extracts of all four tissues of C. didactylus. Five primer combinations were used as described in Katzourakis et al. (2009). The same PCR polymerase employed for mitochondrial bait generation was used, but with an initial denaturing step of 94 °C for 4 min; 35 cycles at 94 °C for 30 s, 55 °C for 30 s, and 72 °C for 30 s; and a final extension of 8 min at 72 °C. An aliquot of each PCR product was visualized on a 1% agarose gel stained with ethidium bromide. To minimize specific biased proviral amplification, PCRs for all five primer combinations were performed in triplicate. PCR products were purified using QIAquick columns (Qiagen) and Sanger sequenced to verify the correct target was amplified. The verified PCR products were quantified using a NanoDrop ND-1000 (Thermo-Scientific) and equimolar products from the five PCRs were pooled. A total of 1.3 µg of the pooled DNA was used as bait for each hybrid capture reaction (Maricic et al. 2010). Additionally, 500 ng of the pool was kept for amplicon Illumina library construction and sequencing

Sample Preparation for Targeted Illumina Sequencing

Aliquots from each DNA extract were used for Illumina library construction. Extinct sloth libraries were constructed in the aDNA facility in a UV laminar hood dedicated to library preparation, while the modern sloth libraries were constructed in a normal DNA laboratory in a different part of the Institute. Additionally, amplicon libraries of SloEFV pol gene amplicons from four tissues of the C. didactylus were constructed together with modern sloth DNA libraries. Library preparation and indexing PCR followed a published protocol (Meyer and Kircher 2010). Each indexed library contained a unique index to allow for subsequent discrimination. The Indexing PCR products were purified using QIAquick columns (Qiagen) with final elution of 40 µl elution buffer (EB).

Triplicate amplifications were performed on each 10 µl eluate of indexing PCR product. PCRs were performed in 50 ml reactions using Herculase II Fusion DNA Polymerase (Agilent Technologies, Catalog 600677) and IS5 and IS6 (Meyer and Kircher 2010) with a final concentration of 400 nM each. Cycling conditions comprised an activation step of lasting 3 min at 95 °C, followed by 20 cycles of denaturation at 95 °C for 20 s, annealing at 60 °C for 25 s, and elongation at 72 °C for 30 s, with a final extension step at 72 °C for 3 min. PCR products were purified using the QIAquick PCR purification kit and measured using Tapstation 2200 (Agilent Technologies Catalog G2964AA). PCR products of the same index were pooled to make sure each had at least 2 µg for hybrid capture following a published protocol (Maricic et al. 2010). The mitochondria and SloEFV pol gene targets were enriched in separate capture reactions for each sample. After 2 days of hybridization and subsequence elution steps, 20 µl EBT was used for final elution using and 1 µl of the final eluate was measured using Tapstation 2200 and 5 µl was taken in triplicate for 18 cycles of amplification using Herculase II Fusion DNA Polymerase following the same protocol as above. The product was purified using the QIAquick PCR purification kit and eluted in 40 µl EB. The eluate was measured using a Tapstation 2200 (Agilent Technologies).

Sequencing of Sloth Mitogenomes and the SloEFV Pol Gene

Four separate pools were generated as follows: The first pool included all mitochondrial enrichment products of modern sloths; the second pool included all mitochondrial enrichment products of extinct sloths; the third pool included all SloEFV pol gene enrichment products of modern sloths; and the fourth pool included all SloEFV pol gene enrichment products of extinct sloths. The four pools were measured on a Tapstation 2200 (Agilent Technologies Catalog G2964AA). Each of the two pools of mitochondrial products was subsequently sequenced on a MiSeq (V2 300 Cycle Kit) at the National High-throughput DNA Sequencing Centre, Copenhagen, Denmark, and each of the two pools of SloEFV gene products was sequenced on a MiSeq (V2 300 Cycle Kit) at the Berlin Center for Genomics in Biodiversity Research, Berlin, Germany.

Sequencing Data Analysis

The sequencing data underwent adaptor trimming with cutadapt-1.2.1 (Martin 2011) using default settings and quality trimming with Trimmomatic-0.22 (Bolger et al. 2014). The minimum size cut-off was set to 20 bp for extinct sloth data and 36 bp for modern sloth data. PCR duplicates (clonality in the sequencing data) with 100% sequence identity were removed using cd-hit-v4.6.1 (Li et al. 2001).

For extant sloth data, mapping was performed using Burrows-Wheeler Aligner (BWA) version 0.7.10 (http://bio-bwa.sourceforge.net/; Li and Durbin 2009) with default setting. Sequencing reads of 2 two-toed sloths were mapped against the C. didactylus (GenBank: AY960980.1) complete mitochondrial genome, and reads of the three-toed sloth were mapped against the B. variegatus mitochondrial genome (NC_006923.1). The 22 published SloEFV pol gene sequences (Katzourakis et al. 2009) and Sanger sequencing results from the pol gene bait were used as references for mapping captured sequence reads. The mapping results were further processed with samtools (Li et al. 2009) and picard (http://broadinstitute.github.io/picard/, last accessed February 23, 2016) for sorting and removal of clonality. The same sequences were used as references for SloEFV pol gene data from extinct sloths.

Because the two extinct sloths in this study are distantly related to any extant sloth (fig. 1), both two- and three-toed sloth mitochondrial genomes were used as references. To achieve the best possible recovery of sequences from the aDNA data, two mapping approaches were used: 1) We followed the recommendations in Schubert et al. (2012) and performed mapping using optimized BWA parameters. The mapping results were visualized in Geneious version 6.18 (http://www.geneious.com/, last accessed February 23, 2016; Kearse et al. 2012). The number of reads that mapped to a reference and length of the reference covered for ancient samples were calculated; 2) We followed a baiting and iterative mapping approach (MITObim version 1.7) using the same reference sequences. For mitochondrial genome reconstruction, the standard procedure (TUTORIAL I) was used, while the quick option (TUTORIAL III) was used for generating SloEFV pol gene consensus sequences. Kbait values for all MITObim analysis on aDNA data were set to 21 bp after we tested a gradient of kbait values for MITObim analysis (supplementary table S1, Supplementary Material online).

Because the MITObim always gives a longer consensus than the BWA approach, we only used MITObim results for further analysis. The two MITObim consensus sequences generated using each of the two modern sloth mitochondrial genome references (C. didactylus and B. variegatus) were aligned using Smith–Waterman criteria (Smith and Waterman 1981) and the resulting consensus of the pairwise alignment was used as reference for remapping the ancient sloth Illumina sequecning reads in Geneious version 6.18. The alignment of the remapping was manually checked by taking 1% of the mapped reads to perform BLASTN in NCBI. The reads matching C. didactylus mitochondrion genome by BLASTN at ≥99% were discarded as contaminants from the bait and all showing ≥99% identity to human mitochondria sequence were discarded as sample contaminants acquired during excavation or later handling. Coverage per base was determined and only positions with at least 8× coverage were retained for further analysis.

Coverage Statistics Plot of Mitochondrial Mapping

Duplicate reads were removed from the mapping file with MarkDuplicates utility from Picard tools (v1.106; http://broadinstitute.github.io/picard/). The coverage of the different genomes was computed with bedtools genomecov (v2.15.0; Quinlan and Hall 2010). The depth of coverage per genome was plotted using R (v3.1.3; R Development Core Team 2008). Each line in the graph is the number of deduplicated reads (per sample) per position in the reference genome.

In Silico Screen for SloEFV Pol Gene Sequences

The assembled C. hoffmanni genome (GenBank ID: ABVD00000000.2) was queried with the Sanger sequencing result of the SloEFV bait amplicon using command line BLASTN (version 2.2.29; Altschul et al. 1990). Genomic sequencing reads from an M. darwinii (SRX327589, SRA, NCBI) sample were downloaded and underwent the same adaptor and quality trimming procedure as described above for aDNA. Mylodon darwinii reads were mapped to the 22 published pol gene partial sequences using MITObim for baiting (fishing for similar reads using a kmer sequence length of 21 bp) and iterative mapping. The resulting consensus sequences from MITObim were used as references for remapping in Geneious version 6.18. The alignments of remapping were checked by eye to ensure there were no gaps and that there was at least a 21 bp overlap among overlapping reads.

Phylogenetic Analyses

We aligned our consensus sequences to annotated mitochondrial genomes of two sloths (C. didactylus and B. variegatus), a vermilinguan (Tamandua tetradactyla), and a cingulate (Dasypus novemcinctus) downloaded from GenBank using the MUSCLE (Edgar 2004) algorithm implemented in Geneious v 8.0.4 (http://www.geneious.com; Kearse et al. 2012). After checking by eye to ensure reading frames were maintained, we deleted the two reference sloth genomes from the alignment to focus analyses on our new sequences, while retaining Dasypus and Tamandua for use as outgroups in subsequent phylogenetic analysis.

We used PartitionFinder v1.1.1 (Lanfear et al. 2012, 2014) to determine the most appropriate partitioning scheme and set of models for our data. We divided the alignment in first, second, and third codon positions for coding regions, 12S, 16S, and tRNAs, but excluded the control region as it was unalignable. Overlapping regions (ATP6/ATP8, ATP6/COIII, ND4/NDL, and ND5/ND6) were assigned to the first gene in the alignment and excluded from the other. We used PartitionFinder’s greedy algorithm, with the Bayesian Information Criterion as our model selection criterion. In considering the balance between model adequacy and overcomplexity, we opted to exclude all models with a proportion of invariant sites (+I) from consideration, as incorporating gamma distributed rates already accounts for this pattern (Allman et al. 2008). The best-fitting scheme resulted in three partitions: 1) Codon position 1 + 12S + 16S + tRNAs; 2) codon position 2; and 3) codon position 3. PartitionFinder assigned a GTR + Γ model of evolution to all partitions except codon position 3, which was assigned an HKY + Γ model.

We performed phylogenetic analyses using Maximum Likelihood (ML) and Bayesian Inference (BI) frameworks. An ML search with 100 rapid bootstrap replicates was performed using RAxML v7.3 (Stamatakis 2006), with the 3 partitions assigned independent GTR + Γ models using the GTRGAMMA option. For analyses under BI, we used MrBayes 3.2.4 (Ronquist et al. 2012), with each partition assigned the best-fit model identified by PartitionFinder. We ran 2 simultaneous analyses, each with 4 chains (1 cold, 3 heated), of 10 million generations, sampling from the chain every 1,000 steps. We ensured that both runs had converged by checking that potential scale reduction factor scores had decreased to ∼1 (Gelman and Rubin 1992), and by visual inspection of trace plots in the program Tracer v 1.6 (http://tree.bio.ed.ac.uk/software/tracer/). After confirming that effective sample sizes were sufficient (>200), we then discarded the first 25% of retained trees as burn-in and summarized the remaining posterior sample using a 50% majority rule consensus tree. Nodal support was summarized as posterior probabilities.

To estimate divergence times of the sampled sloth lineages, we employed a relaxed uncorrelated log-normal molecular clock, as implemented in BEAST v.2.3.1 (Drummond et al. 2012). Because we were primarily interested in inferring divergence times among, and not within, sloth taxa, we replaced our two samples of N. shastensis with a consensus sequence produced in Geneious. Additional analyses including both sequences were performed to ensure that the use of a consensus did not bias divergence time estimates at other nodes in the tree (supplementary table S2, Supplementary Material online). Outgroup taxa were pruned from the alignment for divergence time estimation.

We employed the Fossilized Birth–Death (FBD) process as a tree prior for inferring the timescale of higher level sloth diversification. This tree prior specifically accommodates sampling of noncontemporaneous tips (i.e., fossil taxa) in a robust analytic framework, obviating the need for node age priors that arbitrarily discard fossil observations that are younger than the oldest member of a clade (Gavryushkina et al. 2014; Heath et al 2014). Implementation of the FBD process does not require sampling molecular or morphological data for all extinct taxa. Rather, if we have a sample of fossil taxa whose approximate ages and position on a phylogenetic tree are known, we can integrate over uncertainty regarding their exact topological positions and branching and extinction times while simultaneously sampling parameters of the underlying phylogeny, including divergence times of taxa of interest. As such, the FBD process makes much more complete and efficient use of fossil information than node-based dating for clades such as folivorans that possess relatively good fossil records.

We selected 29 fossil folivorans (supplementary table S3, Supplementary Material online) that could 1) be unambiguously associated with particular higher level groupings on the basis of formal phylogenetic analyses of morphological data and 2) for which relatively well-constrained stratigraphic dates were available. If precise dating is unavailable, we used the rank of the South American Land Mammal Age (SALMA) from which the taxon is known. We supplemented the study with two taxa that have not been subjected to phylogenetic analysis but that are potentially important due to their ages and possible significance. Imagocnus zazae is well constrained temporally to the Burdigalian (16.1–21.5 Ma) of Cuba; if correctly attributed, it is one of the oldest megalonychids recovered to date (MacPhee and Iturralde-Vinent 1995; MacPhee et al. 2003). Pseudoglyptodon chilensis is almost certainly not a crown folivoran, although this can only be stated speculatively as it has not yet been included in a formal analysis. However, the most complete example (and holotype) of this species is temporally constrained to the Tinguirirican SALMA, and more explicitly to 31.5 Ma on the basis of 40Ar/Ar39 and 40K/Ar40 dating (McKenna et al. 2006). This taxon provides useful information about the minimum age of the crown. Full details of fossil taxa and their associated ranges are provided in supplementary tables S3, S4, and S5, Supplementary Material online.

We added these fossil taxa to our molecular matrices, coding each site in the alignment as “?.” We then specified a series of backbone topology constraints (fig. 1) reflecting results from morphology-based phylogenetic studies (Gaudin 1995, 2004; McDonald and Muizon 2002; Rega et al. 2002; Pujos 2006; De Iuliis, Pujos, Bargo, et al. 2009; De Iuliis, Pujos, Tito 2009; De Iuliis et al. 2011; Pujos et al. 2012). These constraints were then enforced during BEAST’s Markov Chain Monte Carlo (MCMC) search. We conditioned our analyses on the proportion of extant sloths sampled in our data which, based on the six xenarthran species (Wilson and Reeder 2005), is 0.5. We further conditioned our analyses on origin time of the FBD process rather than root age, which is more common in node dating analyses. We assigned origin time an offset exponential prior with an offset of 31.5 corresponding to the age of Pseudoglyptodon and a mean = 5, resulting in a 97.5% upper limit of 54.5 Ma. Prior probability distributions are required for three additional parameters under the FBD model: Net diversification rate (speciation minus extinction), turn-over rate (extinction rate divided by speciation rate), and fossil sampling proportion, or the proportion of fossil taxa that are sampled before going extinct (fossil recovery rate divided by the sum of extinction rate and fossil recovery rate). We placed a broad exponential prior with mean = 1 on net diversification rate. For turn-over rate, we used a beta prior with alpha = 2 and beta = 1. This distribution can be thought of as an inverse exponential distribution, and places greatest prior probability on higher turn-over rates (median = 0.7, 95% quantiles = 0.158–0.99), which seems reasonable given the low extant diversity and high fossil diversity of xenarthrans. We also placed a beta prior on fossil sampling proportion, but with alpha = beta = 2. This symmetric prior places greatest weight on intermediate sampling proportions (median = 0.5, 95% quantiles = 0.01–0.9). Two independent Markov Chain Monte Carlo –>analyses were run in BEAST 2.3.1, each of 50 million generations with sampling every 5,000 steps. After ensuring that both chains had converged on the target distribution and that effective sample sizes were greater than 200 for all parameters, we discarded an appropriate number of trees as burn-in (the exact proportion varied among runs but was typically between 10% and 20% ). We combined the remaining samples using LogCombiner v.2.3.1, pruned all fossil taxa except those with molecular data using a custom R script, and produced a maximum clade credibility tree from the remaining sample using TreeAnnotator v.2.3.1.

Two fossil taxa identified as belonging to the folivoran crown in the literature review are of interest. Octodontotherium grande and Deseadognathus riggsi are both known from the Deseadan (29.4–24.2 Ma; Dunn et al. 2012) of Argentina and have been placed in crown families or superfamilies (Mylodontidae/Mylodontonoidea and Megalonychidae, respectively) in formal phylogenetic analyses (Gaudin 1995, 2004; Carlini and Scillato-Yané 2004; Shockey and Anaya 2010). Several factors suggest that these placements should not be accepted unequivocally. First, both taxa are old relative to all other crown folivorans identified in our survey (the next oldest are Eucholoeops and Imagocnus, with maximum ages of 19.33 and 21.5 Ma, respectively MacPhee and Iturralde-Vinent 1995; De Iuliis et al. 2014). Alone, the implication of moderate ghost lineages for crown families is not problematic. However, several authors have argued for a more basal or stem folivoran position for Octodontotherium (alternative hypotheses discussed by Gaudin 2004). Furthermore, the holotype and sole specimen of D. riggsi consists of the rostral portion of a mandible with a caniniform and two molariform teeth, with the result that it could be scored for relatively few characters (Carlini and Scillato-Yané 2004). To explore the impact of alternative phylogenetic hypotheses for the placement of these fossils, we performed three separate divergence time estimations. In the first set of analyses, we followed Gaudin (1995, 2004), Carlini and Scillato-Yané (2004), and Shockey and Anaya (2010) in treating the Deseadan fossils as crown taxa, constraining them to fall within Mylodontidae and Megalonychidae, respectively, and thus ensuring a minimum Late Oligocene age for crown folivorans. In the second analysis, we treated the Deseadan taxa as stem folivorans by constraining them to fall outside of crown Folivora. In the third and final analysis, we allowed for the possibility that Deseadan taxa belong to the crown group but not to any of the recognized crown families by allowing them to float between crown and stem positions. Median node ages and associated posterior distributions for each set of analyses were compared to obtain a spectrum of possible divergence times based on alternative phylogenetic histories for the earliest fossil folivorans.

Phylogenetic inference from aligned SloEFV pol gene sequences was conducted in RAxML v8.2.0, using ten independent starting points and retaining the best tree. We evaluated node support from 100 rapid bootstrap replicates.

Sequence Accession

Mitochondrial DNA consensus sequences generated in this study are accessible from GenBank for C. didactylus, C. hoffmanni, B. variegatus, M. darwinii, and N. shastensis (accession numbers KR336791, KR336792, KR336793, KR336794, and KR336795, respectively). SloEFV consensus sequences generated in this study can be accessed in GenBank (accession numbers KR703280–KR703472). Illumina sequence data from which the consensus sequences are derived were deposited in the Sequence Read Archive (SRA) (SRR2007671–SRR2007675).

Results

Hybridization Capture Sequence Retrieval from Extinct and Extant Sloths

MITObim is a new algorithm which iteratively allows for the building of mitochondrial genomes without a closely related or full reference sequence (Hahn et al. 2013). This is particularly crucial for analysis of the extinct sloths: The extant sloths are quite distantly related to the two extinct sloth clades in this study, raising the possibility that mapping results using BWA only could be biased or incomplete, resulting in loss of data. Comparing mapping results from the fossil samples using MITObim and BWA, MITObim achieved a 10- to 15-fold increase in the number of mappable reads and a subsequent 12- to 17-fold increase in total mitogenome coverage (table 2). After preprocessing, the three modern sloths were mapped using BWA to either a B. variegatus or C. didactylus reference genome, yielding as expected, much higher coverage and enrichment efficiency than ancient samples (table 2).

Table 2

Mapping Statistics

Sample	Total Sequences	Sequences Mapped with BWA	Sequences Mapped with MITObim	Percentage of Total Sequences Mapped with BWA	Percentage of Total Sequences Mapped with MITObim	BWA Contig Length (bp)	MITObim Contig Length (bp)	Percentage of Reference Covered by BWA Mapping	Percentage of Reference Covered by MITObim Mapping
Mylodon darwinii (bone)	2,254,520	617	6,970	0.03	0.31	1,153	15,792	6.97	95.46
M. darwinii (coprolite)	93,426	4	12	0.00	0.01	32	154	0.19	0.93
Nothrotheriops shastensis (Gypsum cave)	912,866	831	47,733	0.09	5.23	997	14,205	6.03	85.87
N. shastensis (Rampart cave)	1,101,424	92	679	0.01	0.06	269	3,796	1.63	22.95
Bradypus variegatus	2,560,964	1,115,396		43.55%
Choloepus didactylus	2,661,464	1,573,533		59.12%
Choloepus hoffmanni	3,796,487	374,716		9.87%
Negative control	6,572	2	8	0.03	0.12	27	37	0.16	0.22

Note.—The C. didactylus mitochondrion genome (NC_006924.1) was used as a mapping reference.

SloEFV Pol Gene Retrieval

A total of 198 SloEFV pol gene sequences were obtained from extant and extinct sloths. Hybridization capture yielded 21 pol genes from B. variegatus, 18 pol genes from C. hoffmanni, and 109 from the 2 C. didactylus (one used for mitochondrial sequencing, the other used for generating bait), while among the extinct sloths, SloEFV sequences could only be retrieved from M. darwinii, and 12 sequences fully covering the reference were identified using MITObim. In silico search yielded 22 additional SloEFV pol gene sequences from published C. hoffmanni genome. Additionally, 16 SloEFV pol gene sequences were retrieved from the Mylodon SRA.

Sloth Mitochondrial Phylogeny

Analyses under ML and BI resulted in identical topologies that confirm results from previous aDNA studies. We recovered two reciprocally monophyletic sloth clades (fig. 2), one containing M. darwinii and the two Choloepus species, the other B. variegatus and N. shastensis. All nodes were unambiguously supported (bootstrap = 100%, PP = 1.0).

F

Time-calibrated phylogeny of extant and extinct sloths inferred from mitogenomic data. Divergence times illustrated represent median ages, while the light blue bars give the 95% HPD intervals for each node. Solid lines give the divergence times inferred when treating Deseadan fossil taxa as representatives of crown families. Dashed lines indicate how median ages change when forcing Deseadan taxa to the folivoran stem (see text for details).

Divergence time estimates are, as expected, highly sensitive to the treatment of Deseadan fossil taxa (table 3). If these fossils are constrained to fall within crown families, as suggested by phylogenetic analysis of morphological data, then we automatically recover a Late Oligocene (∼29 Ma) crown age for Folivora. However, treating these fossils as either stem taxa or allowing them to be crown taxa that fall outside of recognized families reduces the age of crown Folivora by ∼8.5 Ma, bringing it into the Early Miocene, albeit with a lower 95% Highest posterior density interval (HPD) > bound that includes the latest Oligocene. Inter-familial divergences are similarly affected. The divergence between Megalonychidae and Mylodontidae is, as expected, estimated to have occurred in the Late Oligocene (∼27.5 Ma), when the Deseadan taxa are treated as belonging to these families, but reduces to the late Early Miocene (∼19 Ma) when using alternative treatments. The divergence of Bradypus from Nothrotheriidae is estimated to be latest Oligocene/earliest Miocene (23.5 Ma), with a 95% HPD interval spanning Late Oligocene to Middle Miocene, when treating Deseadan fossils as belonging to crown families, but is restricted to the Early Miocene (∼17.5 Ma) when using alternative placements. The divergence of the two Choloepus species is least affected by alternative phylogenetic treatments of the Deseadan fossils, with median estimates of 7.5–5.6 Ma but wide HPD intervals in all cases that span Late Miocene to Early Pleistocene (table 3).

Table 3

Median (95% highest posterior density interval) Age Estimates for Divergence Events in Sloth Phylogeny based on Three Alternative Treatments of Deseaden Taxa

Node	Deseaden Crown	Deseaden Stem	Free
Crown Folivora	29.3 (24.5–34.6)	20.6 (18.0–24.3)	20.8 (17.3–25.7)
Bradypus + Nothrotheriops	23.6 (15.22–31.5)	17.7 (12.9–23.1)	17.3 (12.3–23.4)
Choloepus + Mylodon	27.5 (23.7–31.5)	19.2 (17.6–22.4)	19.5 (16.8–23.3)
Choloepus didactylus + Choloepus hoffmanni	7.5 (3.1–12.6)	5.6 (2.5–9.0)	5.4 (1.4–8.6)

Comparison of SloEFV Phylogeny to Sloth Phylogeny

The topology of the tree based on SloEFV pol gene sequences from extant and extinct sloth species (fig. 3) is not congruent with previous efforts to derive sloth phylogeny from mitochondrial evidence, but it is in good agreement with the tree inferred for modern sloths by Katzourakis et al. (2009). For example, most clades were either exclusively Choloepus derived or grouped SloEFVs from Choloepus and Bradypus (fig. 3). Comparison of the expanded SloEFV phylogeny with reference to the well-resolved sloth host phylogeny clearly demonstrated multiple independent introgression and expansion events at distinct points in sloth evolutionary history. Lineages were found in Choloepus and Bradypus, but not Mylodon (sequences of which are highlighted in red and brown). Absence from Mylodon could represent either a methodological issue (e.g., the poor quality of the sample not allowing for retrieval of this clade) or species-specific resistance to or loss of SloEFV as has been observed even among closely related taxa (e.g., the Pan troglodytes ERV [PtERV] which is common in chimpanzees but absent from humans; Mun et al. 2014).

F

Phylogeny of SloEFV partial pol genes from extant and extinct sloth. The phylogenetic tree was calculated using RAxML; sequences from different sloth species are designated by color. Bootstrap support is indicated by circles 75–90 = white; 90–99 = gray; 100 = black. Sequences for pol genes were obtained for Bradypus variegatus (dark orange), Choloepus hoffmanni (dark blue), Choloepus didactylus (aquamarine), and Mylodon darwinii (red) by hybridization capture. Sequences for B. variegatus and Bradypus pygmaeus (light orange), C. hoffmanni (light blue) and M. darwinii (brown) were obtained in silico. Denoted in yellow are the pol gene sequences used as outgroups for tree construction: Bovine (NC_001831.1), Equine (NC_002201.1), Feline (NC_001871.1), Human (GenBank: Y07725.1), and Rhinolophus (GenBank: JQ814855.1) foamy viruses (all downloaded from NCBI).

The two lineages highlighted in light red, light blue, and aquamarine or brown and light blue represent SloEFVs common to Choloepus and Mylodon only. They represent SloEFV invasion events occurring after the Choloepus–Mylodon clade diverged from other sloth lineages but predate divergence between the families Megalonychidae and Mylodontidae. The lineages highlighted in any combination of dark blue, light blue, and aquamarine were only identified within genus Choloepus and represent relatively recent SloEFV invasions and expansion.

Discussion

Performance of Hybridization Capture

In our experiments, hybridization capture achieved higher enrichment efficiency for modern than ancient samples (table 3), which is expected given the highly fragmented nature of aDNA, contamination with DNA from microorganisms post mortem and the low DNA concentrations in the samples. Nonetheless, hybridization capture produced a roughly 16-fold increase in DNA sequence information for the two extinct sloths compared with published genetic information obtained by PCR (Greenwood et al. 2000), which yielded relatively short sequences. Similar enrichment was observed for SloEFV sequences. Although the enrichment efficiency, SloEFV sequences retrieved compared with total sequences, of the retroviral gene from Mylodon was only 0.002%, the present effort marks the first successful retrieval of retroviral sequences from a Pleistocene sample using hybridization capture.

MITObim provided an additional benefit to sequence retrieval analysis. BWA is a standard mapping software which performs well for mapping low-divergent sequences against a reference genome of the target species or a very closely related species. In contrast, MITObim is able to reconstruct genomes of limited size, like those characteristic of mitochondria, by using either short seed sequences developed from the target species itself or a reference gene/genome from a more distantly related species as a starting reference. Given the deep divergence estimates for sloth clades, it was anticipated that a significant amount of sequence information might be lost if BWA were used, because of both the short length of aDNA sequences and substantial sequence divergence among and between extinct and extant sloth lineages. MITObim increased sequence retrieval from less than 10% to 95.46% and 85.87% for M. darwinii and N. shastensis, respectively. MITObim was also effective in fully covering the targeted SloEFV pol sequence, indicating that the algorithm is useful for accessing nonmitochondrial loci in general.

Toward a Compatible Molecular/Morphological Phylogeny of Folivora

Although studies of the molecular phylogenetics of major mammalian taxa are now commonplace (Murphy et al. 2001; Springer et al. 2011; Dos Reis et al. 2012), some groups remain poorly investigated. Among these are the sloths, the folivoran xenarthrans. As recently as the end of the Pleistocene, this clade included a rich diversity distributed through much of South America, parts of the West Indies, and southern North America. Of this widespread fauna only two genera remain, the tree sloths Bradypus (three-toed sloth) and Choloepus (two-toed sloth); the relationships of these unlikely survivors to their manifold precursors have long engendered controversy in both morphological and molecular studies (figs. 1 and 2).

In this study, we recovered a well-supported sister group relationship between the three-toed sloth (B. variegatus) and the Pleistocene nothrotheriid ground sloth N. shastensis, consistent with previous molecular studies (Greenwood et al. 2000; Poinar et al. 2003). Although recent phylogenetic analyses of morphological character data have tended to place Bradypus as the sister lineage to all other living and extinct sloths (Gaudin 1995, 2004; but see White and MacPhee 2001), the view that Bradypus is related to megatheriid and nothrotheriid ground sloths in general, and to nothrotheriid ground sloths in particular, has a long history in the paleontological literature (Guth 1961; Patterson and Pascual 1968; Webb 1985; Patterson et al. 1992). The precise placement of Bradypus with respect to megatheriids and nothrotheriids remains enigmatic, however, because many apparently autapomorphic character states in Bradypus are also convergently derived in Megatheriidae (Planops, Megatherium, and Eremotherium; Gaudin 1995, 2004). Alternatively, Bradypus may ultimately turn out to be the sister lineage to all other megatheriid and nothrotheriid sloths. At present, Nothrotheriops is the only fossil megathere for which aDNA data are available, although Buckley et al. (2015) have recently published collagen proteomic sequence data that suggest that Megatherium is not a close relative of Bradypus. In that study the two surviving genera of sloths appear as each other’s closest relative—a notable departure from the majority of recent studies cited above. Although our analyses are consistent with the idea that the smallest extant sloths (B. pygmaeus, 2.2–3.5 kg; Anderson and Handley 2002) form a clade with the largest sloths to have ever lived (Megatherium americanum, >3,500 kg; Fariña et al. 1998), full resolution of relationships must therefore await more data.

Choloepus presents other puzzles. Since the mid-20th century most workers have placed Choloepus in Megalonychidae, and a large number of morphological characters, including features of the ear region, palate, snout, and, most strikingly, the caniniform first upper and lower teeth, appear to support this placement (Patterson and Pascual 1968; Webb 1985; Patterson et al. 1992; Gaudin 1995, 2004; White and MacPhee 2001; Pujos et al. 2007, 2012, McDonald et al. 2013). The relationship of Megalonychidae to other sloth clades has proved more difficult to resolve. Based on an analysis of the ear region, Gaudin (1995) recovered a clade comprised of Mylodontidae and Megalonychidae, but after adding cranial and dental characters he (Gaudin 2004) suggested a sister relationship between Megalonychidae and Megatheriidae + Nothrotheriidae. Molecular data are consistent with the signal from the ear region in supporting a sister relationship between Megalonychidae and Mylodontidae (Höss et al. 1996; Greenwood et al. 2000). The conflicting phylogenetic signal found among morphological partitions is particularly concerning because many of the earliest putative representatives of the three crown sloth families are represented by fragmentary material. As already noted, the earliest putative megalonychid D. riggsi from the Deseadan of Patagonia is diagnosed primarily on tooth shape and the size of the diastema between the caniniform and first molariform (Carlini and Scillato-Yané 2004). If some of these features represent symplesiomorphies that are retained in individual crown families, then our ability to resolve phylogenetic relationships and estimate divergence times (see below) at the base of the crown folivoran radiation may be severely compromised. Results of our molecular phylogenetic analysis provide additional support for the emerging view of a mylodontid–megalonychid clade and will, we hope, help future morphological investigations aimed at discovering characters beyond the ear region that support this relationship.

Our molecular estimates for the timescale of higher level folivoran diversification proved highly sensitive to the placement of Deseadan fossil taxa. This in turn yields dramatically different interpretation regarding the ecological context in which sloth diversification occurred. Treating Deseadan taxa as early representatives of crown families, as morphological phylogenetic studies have suggested (Gaudin 1995, 2004; Carlini and Scillato-Yané 2004), generates an Oligocene estimate for the origin of crown Folivora, with interfamilial diversification occurring in the Late Oligocene or earliest Miocene. Most sloth families have their first appearance in Santacrucian or younger deposits (Carlini and Scillato-Yané 2004; Pujos and De Iuliis 2007; Bargo et al. 2012), the so-called “heyday of nonmylodontid diversity” (Webb 1985). As such, an Oligocene origin would be suggestive of a long-fuse model for sloth diversification, with initial family-level divergences occurring some 5–10 Myr prior to the dramatic generic diversification indicated in the fossil record. Alternatively, treating Deseadan fossils as stem folivorans generates an origin of crown Folivorans at the Oligocene/Miocene boundary, with interfamilial diversification occurring rapidly in the Early Miocene. Delsuc et al. (2004) recovered a similar Early Miocene estimate for the divergence of extant sloth genera and, noting the correspondence between this date and the proposed radiation of tolypeutine armadillos, suggested a role for the onset of rainshadow conditions east of the Andes and a transition to dryer, more open environments in initiating the crown sloth radiation (Blisniuk et al. 2005).

Despite uncertainty regarding the exact timescale of folivoran diversification, relative ages suggest that the origins and higher level diversification of crown group sloths occurred within a narrow temporal window. The conspicuous lack of consensus in efforts to resolve higher level folivoran relationships may therefore in part be the result of extremely rapid diversification of major sloth lineages that, in turn, provided little time for the accumulation of multiple unambiguous morphological synapomorphies. Our analyses suggest that crown folivorans originated and rapidly gave rise to at least five families of Bradypodidae, Megatheriidae, Nothrotheriidae, Mylodontidae, and Megalonychidae within ∼3–6 Myr during the Late Oligocene/Early Miocene. If the Deseadan taxa used in our computations are in fact stem folivorans, they might have exhibited mosaics of symplesiomorphies that were subsequently retained in individual crown clades, making it hard to draw a line between them. The correct placement of Orophodontinae (or Orophodontidae) is similarly at issue, with several authors having suggested a more stem-ward position for this taxon outside of Mylodontidae or crown group Folivora (discussed by Gaudin 2004). Future phylogenetic work and more precise age estimates may ultimately pull the age of crown group Folivora back into the Oligocene but, as the Early Miocene Santacrucian record attests, the rapid diversification of crown families is likely a real feature of the data and will presumably present challenges to resolving the earliest stages of this radiation. Just as detailed examination of the morphology of Early Miocene taxa has helped to elucidate their phylogenetic relationships (e.g., Eucholoeops; Gaudin et al. 2015), reexamination of Oligocene taxa for which representative material is available, such as Octodontotherium, should be a priority for xenarthran systematic studies in the coming years.

Perhaps the most surprising result to emerge from our analyses is the finding that the two extant species of Choloepus diverged ∼6–7 Ma. Previous attempts to understand the timing of xenarthran diversification have focused on higher level divergences (Sarich 1985; Delsuc et al. 2004) and few molecular data have been available for C. hoffmanni until now. Of special interest in this regard is the recent estimate of 12 Ma for the divergence of B. torquatus from B. variegatus and B. tridactylus (Moraes-Barros et al. 2011), which indicates that an ancient divergence time for the two Choloepus species should not be unexpected. There are also a number of striking morphological differences between C. hoffmanni and C. didactylus, including the number of cervical and caudal vertebrate (Wetzel 1985; Maslin et al. 2007; Hautier et al. 2010; Asher et al. 2011). However, our divergence time estimates are particularly remarkable given that Choloepus itself has recently been nested within the endemic Antillean megalonychid radiation as a close relative of Acratocnus (Gaudin 1995, 2011, 2004; White and MacPhee 2001, McDonald et al. 2013, but see Pujos et al. 2007; Lyon et al. 2015). Although the earliest firm evidence for an Antillean megalonychid is I. zazae from the Early Miocene (Burdigalian, 16.1–21.5 Ma) site of Domo de Zaza, Cuba (MacPhee et al. 2003), virtually all other West Indian fossil material comes from Late Pleistocene and Holocene sites, with the result that the timescale of sloth diversification among and within these islands is currently very poorly understood. It is conceivable, if unlikely, that ancestral Choloepus (re)colonized the mainland from the islands in the Late Miocene. Alternatively, Choloepus may be more distantly related to the Antillean radiation (Pujos et al. 2007). An intraspecific analysis based on the mitochondrial gene cytochrome oxidase I lends some support to this idea; Moraes-Barros and Arteaga (2015) estimate a similarly ancient divergence (∼7 Ma, based on a strict molecular clock) between Central and South American populations of C. hoffmanni. This finding suggests that Choloepus’s disjunct present day distribution was affected by uplift of the northern Andes and supports the hypothesis of a long history on the mainland.

SloEFV Macroevolution with Sloths

ERVs are expected to follow a pattern of virus–host coevolution. Once endogenized and fixed in a species, retroviruses adopt the relatively slow mutation rate of the host and from that point onward, host and viral phylogenetic relationships should correlate. However, subsequent viral invasion by related exogenous retrovirus or intracellular retrotransposition and expansion can lead to a pattern of evolution that is highly discordant with host phylogeny and more closely resembles waves of infection (Hayward et al. 2013, 2015). A mixture of both of these evolutionary trajectories can be observed in the reconstructed history of SloEFV endogenization in sloths. Sequences were observed that predate the predicted divergence of any of the sloth lineages examined in this study. These occupy basal positions in the SloEFV tree (fig. 3) and must reflect ancient retroviral colonization events occurring in the last 100 Myr (Katzourakis et al. 2014). However, many of sequences obtained in this study were lineage specific, with Choloepus exhibiting an abundance of private lineages not found in other extant or extinct sloth groups. There were also several cases of SloEFV sequences demonstrating close affinity of Choloepus and Bradypus, in contrast with the mitochondrial analysis which clearly demonstrated a close phylogenetic relationship between Choloepus and Mylodon (figs. 2 and 3). There are several likely reasons for discordance. The degraded nature of the Mylodon DNA makes confirmation of the absence of a given sequence group problematic. For example, no SloEFV sequences could be obtained from Nothrotheriops although mitochondrial DNA sequences were obtained. Similarly, SloEFV lineages in the two Mylodon samples in the current study differed to some extent, suggesting incomplete sample coverage. Another explanation could be variation in resistance to the exogenous ancestor of the SloEFV, resulting in the presence of the observed lineages in some taxa but absence in others. The chimpanzee ERV PtERV is common in the genomes of gorillas and chimpanzees but absent in humans, despite overwhelming molecular support for humans and chimpanzees being each other’s closest relative, to the exclusion of gorillas (Mun et al. 2014). Similar resistance to specific SloEFVs could have prevented acquisition of specific SloEFVs in Mylodon. Finally, polymorphism in repetitive elements, such as Alu, can lead to the loss of specific retroelement lineages in different taxa, thus mimicking lineage-specific viral resistance (Salem et al. 2003).

In light of our divergence estimates for different sloth clades, the abundance of lineages specific to Choloepus–Bradypus, Mylodon–Choloepus, and Choloepus–Choloepus implies that SloEFVs have been actively spreading and amplifying in sloth lineages during the last 20 Myr, with activity in the 10 Myr since the divergence of the Choloepus/Mylodon lineages being particularly marked. Future comparisons of the full genomes of extinct and extant sloths should focus on repetitive elements to discern whether relative permissiveness in retroviral lineage accumulation is a common feature of retroelements in Choloepus or is specific to SloEFV.

Supplementary Material

Supplementary tables S1–S5 and references are available at Genome Biology and Evolution online (http://www.gbe.oxfordjournals.org/).