Maxime Courcelle
email: maxime.courcelle@umontpellier.fr

Marie-Ka Tilak

Yuri L R Leite

Emmanuel J P Douzery

Pierre-Henri Fabre

Digging for the spiny rat and hutia phylogeny using a gene capture approach, with the description of a new mammal subfamily

Keywords: Rodents, Echimyidae, Capromyinae, Exon capture, Phylogenomics

Next generation sequencing (NGS) and genomic database mining allow biologists to gather and select large molecular datasets well suited to address phylogenomics and molecular evolution questions. Here we applied this approach to a mammal family, the Echimyidae, for which generic relationships have been difficult to recover and often referred to as a star phylogeny. These South-American spiny rats represent a family of caviomorph rodents exhibiting a striking diversity of species and life history traits. Using a NGS exon capture protocol, we isolated and sequenced ca. 500 nuclear DNA exons for 35 species belonging to all major echimyid and capromyid clades.

Exons were carefully selected to encompass as much diversity as possible in terms of rate of evolution, heterogeneity in the distribution of site-variation and nucleotide composition.

Supermatrix inferences and coalescence-based approaches were subsequently applied to infer this family's phylogeny. The inferred topologies were the same for both approaches, and support was maximal for each node, entirely resolving the ambiguous relationships of previous analyses. Fastevolving nuclear exons tended to yield more reliable phylogenies, as slower-evolving sequences were not informative enough to disentangle the short branches of the Echimyidae radiation. Based on this resolved phylogeny and on molecular and morphological evidence, we confirm the rank of the Caribbean hutias -formerly placed in the Capromyidae family -as Capromyinae, a clade nested within Echimyidae. We also name and define Carterodontinae, a new subfamily of Echimyidae, comprising the extant monotypic genus Carterodon from Brazil, which is the closest living relative of West Indies Capromyinae.

INTRODUCTION

Over the past years, the advent of next-generation DNA sequencing (NGS) shifted the way biologists access and analyze molecular data [START_REF] Zhang | The impact of next-generation sequencing on genomics[END_REF][START_REF] Goodwin | Coming of age: ten years of next-generation sequencing technologies[END_REF]. Much emphasis has been put on the vast amount of data NGS can provide in a broad range of organisms.

NGS are now routinely used in diverse biological fields for metabarcoding [START_REF] Bohmann | Environmental DNA for wildlife biology and biodiversity monitoring[END_REF], population genomics [START_REF] Yi | Sequencing of 50 Human Exomes Reveals Adaptation to High Altitude[END_REF], ancient DNA studies [START_REF] Rohland | Proboscidean Mitogenomics: Chronology and Mode of Elephant Evolution Using Mastodon as Outgroup[END_REF]) and phylogeny [START_REF] Misof | Phylogenomics resolves the timing and pattern of insect evolution[END_REF] purposes. Thus, NGS offers the possibility to choose molecular markers with desirable properties, rather than being limited by commonly used markers. This approach allows biologists to use the most appropriate molecular markers for their research, be it single nucleotide polymorphism, ultra-conserved elements, or genes [START_REF] Peterson | Double Digest RADseq: An Inexpensive Method for De Novo SNP Discovery and Genotyping in Model and Non-Model Species[END_REF][START_REF] Faircloth | Ultraconserved Elements Anchor Thousands of Genetic Markers Spanning Multiple Evolutionary Timescales[END_REF][START_REF] Yang | Clinical Whole-Exome Sequencing for the Diagnosis of Mendelian Disorders[END_REF]). Among NGS-based techniques, sequence capture uses baits built from genomic references to specifically enrich sequencing libraries in genomic regions of interest. In the field of phylogenomics, sequence capture has been used to reliably gather datasets of hundreds to thousands of orthologous targets among various clades, such as mammals [START_REF] Schweizer | Targeted capture and resequencing of 1040 genes reveal environmentally driven functional variation in grey wolves[END_REF], birds [START_REF] Smith | Target Capture and Massively Parallel Sequencing of Ultraconserved Elements for Comparative Studies at Shallow Evolutionary Time Scales[END_REF], squamates [START_REF] Bragg | Phylogenomics of a rapid radiation: the Australian rainbow skinks[END_REF], lissamphibians [START_REF] Mccartney-Melstad | Exon capture optimization in amphibians with large genomes[END_REF][START_REF] Newman | Sequence capture and next-generation sequencing of ultraconserved elements in a large-genome salamander[END_REF] and gymnosperms [START_REF] Suren | Exome capture from the spruce and pine giga-genomes[END_REF]. Several methodological benefits of sequence capture can explain this recent interest in phylogenomics. Firstly, datasets as large as whole exomes can be acquired relatively cheaply [START_REF] Bi | Transcriptome-based exon capture enables highly cost-effective comparative genomic data collection at moderate evolutionary scales[END_REF]. Secondly, capture performs well on degraded material, which enable the use of DNA from museum specimens [START_REF] Fabre | Rodents of the Caribbean: origin and diversification of hutias unravelled by nextgeneration museomics[END_REF][START_REF] Mccormack | Sequence capture of ultraconserved elements from bird museum specimens[END_REF]. Thirdly, the technique tolerates a certain amount of nucleotide variation between baits and captured sequences, usually 15% of molecular divergence between a target and its reference [START_REF] Bi | Transcriptome-based exon capture enables highly cost-effective comparative genomic data collection at moderate evolutionary scales[END_REF]. Thus, sequences from non-model or even extinct species might be captured if close relatives are available in genomic databases [START_REF] Christmas | Finding needles in a genomic haystack: targeted capture identifies clear signatures of selection in a nonmodel plant species[END_REF][START_REF] Carpenter | Pulling out the 1%: Whole-Genome Capture for the Targeted Enrichment of Ancient DNA Sequencing Libraries[END_REF][START_REF] Enk | Ancient Whole Genome Enrichment Using Baits Built from Modern DNA[END_REF].

Sequence capture offers the possibility to choose a dataset that fits the question being studied. Indeed, several characteristics linked to molecular evolution parameters have indeed proven useful to resolve difficult nodes of the phylogeny of metazoan organisms [START_REF] Rodríguez-Ezpeleta | Detecting and Overcoming Systematic Errors in Genome-Scale Phylogenies[END_REF][START_REF] Philippe | Resolving Difficult Phylogenetic Questions: Why More Sequences Are Not Enough[END_REF][START_REF] Romiguier | Less Is More in Mammalian Phylogenomics: AT-Rich Genes Minimize Tree Conflicts and Unravel the Root of Placental Mammals[END_REF]. Firstly, the evolutionary rate of molecular markers is defined by the number of substitutions that DNA, RNA and protein sequences accumulates through time. It thus impacts saturation and more generally homoplasy rates [START_REF] Philippe | How good are deep phylogenetic trees?[END_REF]. As such, a number of authors recommended to only use slower-evolving markers in phylogenetic studies [START_REF] Steppan | Phylogeny and Divergence-Date Estimates of Rapid Radiations in Muroid Rodents Based on Multiple Nuclear Genes[END_REF][START_REF] Betancur-R R | Addressing gene tree discordance and non-stationarity to resolve a multi-locus phylogeny of the flatfishes (Teleostei: Pleuronectiformes)[END_REF][START_REF] Nosenko | Deep metazoan phylogeny: When different genes tell different stories[END_REF].

However, such slower-evolving genes may not accumulate enough substitutions to be phylogenetically informative [START_REF] Townsend | Phylogenetic signal and noise: predicting the power of a data set to resolve phylogeny[END_REF]. Substitutions in conserved genes may also be under selective constraints themselves and, as a result, be more prone to homoplasy [START_REF] Källersjö | Homoplasy Increases Phylogenetic Structure[END_REF]). This led many recent studies to combine both faster-and slower-evolving genes in order to benefit from their resolution at different taxonomic scales (de Knijff 2000; [START_REF] Delsuc | Resolving the phylogenetic position of Darwin's extinct ground sloth (Mylodon darwinii) using mitogenomic and nuclear exon data[END_REF].

Secondly, distribution of substitutions across a gene may also have an impact on the number of homoplasies. If a site is more frequently subject to multiple hits, it is more likely to reach saturation, and to undergo convergences and/or reversals. Ideally, the variability of a phylogenetic marker should then be distributed equally among sites [START_REF] Townsend | Profiling Phylogenetic Informativeness[END_REF]. An example of such a marker is the BReast CAncer 1 (BRCA1) gene, for which all codon positions were shown to accumulate substitutions at a similar rate and thus was a useful marker in mammal phylogeny [START_REF] Delsuc | Molecular Phylogeny of Living Xenarthrans and the Impact of Character and Taxon Sampling on the Placental Tree Rooting[END_REF][START_REF] Meredith | Impacts of the Cretaceous Terrestrial Revolution and KPg Extinction on Mammal Diversification[END_REF]. Thirdly, the nucleotide composition of sequences has been linked to their potential to infer the correct phylogeny [START_REF] Romiguier | Analytical Biases Associated with GC-Content in Molecular Evolution[END_REF]. G+C-rich regions have for example been correlated to recombination through biased gene conversion [START_REF] Galtier | GC-content evolution in mammalian genomes: the biased gene conversion hypothesis[END_REF]. By reducing the available substitution pool, biased gene conversion mechanically increases the rate of homoplasy [START_REF] Romiguier | Analytical Biases Associated with GC-Content in Molecular Evolution[END_REF]. Moreover, recombination hotspots are regions where ancestral polymorphism is more often artifactually sorted, which leads to conflicts between gene trees and species tree [START_REF] Hobolth | Incomplete lineage sorting patterns among human, chimpanzee, and orangutan suggest recent orangutan speciation and widespread selection[END_REF]. Some authors thus recommend favoring A+T-rich genes as phylogenetic markers [START_REF] Betancur-R R | Addressing gene tree discordance and non-stationarity to resolve a multi-locus phylogeny of the flatfishes (Teleostei: Pleuronectiformes)[END_REF][START_REF] Romiguier | Less Is More in Mammalian Phylogenomics: AT-Rich Genes Minimize Tree Conflicts and Unravel the Root of Placental Mammals[END_REF]).

As such, the lack of diversity of markers in phylogenetic analyses may hamper or at least affect inferences.

In recent phylogenomic studies, sequence capture has primarily targeted ultra-conserved elements (UCE) [START_REF] Smith | Target Capture and Massively Parallel Sequencing of Ultraconserved Elements for Comparative Studies at Shallow Evolutionary Time Scales[END_REF][START_REF] Manthey | Comparison of Target-Capture and Restriction-Site Associated DNA Sequencing for Phylogenomics: A Test in Cardinalid Tanagers (Aves, Genus: Piranga)[END_REF][START_REF] Mccormack | Sequence capture of ultraconserved elements from bird museum specimens[END_REF][START_REF] Newman | Sequence capture and next-generation sequencing of ultraconserved elements in a large-genome salamander[END_REF][START_REF] Branstetter | Enriching the ant tree of life: enhanced UCE bait set for genome-scale phylogenetics of ants and other Hymenoptera[END_REF]). These slow-evolving regions are mostly conserved between species, making them an ideal target for sequence-capture baits. The phylogenetic signal carried by UCEs and their flanking regions has also been shown to be informative [START_REF] Giarla | The Challenges of Resolving a Rapid, Recent Radiation: Empirical and Simulated Phylogenomics of Philippine Shrews[END_REF][START_REF] Gilbert | Genomewide ultraconserved elements exhibit higher phylogenetic informativeness than traditional gene markers in percomorph fishes[END_REF]. However, studies at lower taxonomic levels depend almost entirely on UCE flanking regions in order to gather enough variability between closely related species [START_REF] Smith | Target Capture and Massively Parallel Sequencing of Ultraconserved Elements for Comparative Studies at Shallow Evolutionary Time Scales[END_REF]. Unlike UCEs [START_REF] Bejerano | Ultraconserved Elements in the Human Genome[END_REF][START_REF] Katzman | Human Genome Ultraconserved Elements Are Ultraselected[END_REF], nuclear exons cover a wider range of molecular evolution properties [START_REF] Ranwez | OrthoMaM: A database of orthologous genomic markers for placental mammal phylogenetics[END_REF], while being conserved enough for a sequence capture kit to be designed [START_REF] Bragg | Exon capture phylogenomics: efficacy across scales of divergence[END_REF][START_REF] Hugall | An Exon-Capture System for the Entire Class Ophiuroidea[END_REF]. Thus, nuclear exons seem well-suited molecular markers to resolve difficult phylogenetic questions.

The generic relationships among the South-American spiny rats (Echimyidae), as well as their affinities with Caribbean hutias (Capromyidae) present such a challenging phylogenetic problem. Indeed, the phylogenetic relationship of hutias from the Greater Antilles have been widely debated and remain a major topic in rodent systematics [START_REF] Woods | Adaptive Radiation of Capromyid Rodents: Anatomy of the Masticatory Apparatus[END_REF][START_REF] Mckenna | Classification of mammals: above the species level[END_REF][START_REF] Woods | Insular patterns and radiations of West Indian rodents[END_REF][START_REF] Fabre | Rodents of the Caribbean: origin and diversification of hutias unravelled by nextgeneration museomics[END_REF]2017). Earliest classifications of these rodent families, using molecular or morphological data, were incongruent and mostly highlighted their very fast early diversification [START_REF] Lara | The Simultaneous Diversification of South American Echimyid Rodents (Hystricognathi) Based on Complete Cytochrome b Sequences[END_REF][START_REF] Leite | Evolution of South American spiny rats (Rodentia, Echimyidae): the star-phylogeny hypothesis revisited[END_REF][START_REF] Carvalho | Relationships among extant and fossil echimyids (Rodentia: Hystricognathi)[END_REF][START_REF] Emmons | A Revision of the Genera of Arboreal Echimyidae (Rodentia: Echimyidae, Echimyinae), With Descriptions of Two New Genera[END_REF]. Capromyidae have then been classified either as a family or as a subfamily, considering both morphological [START_REF] Woods | Adaptive Radiation of Capromyid Rodents: Anatomy of the Masticatory Apparatus[END_REF][START_REF] Mckenna | Classification of mammals: above the species level[END_REF][START_REF] Verzi | Phylogeny and Evolutionary Patterns of South American Octodontoid Rodents[END_REF]) and molecular inferences [START_REF] Upham | Evolution of caviomorph rodents: a complete phylogeny and timetree for living genera[END_REF][START_REF] Fabre | Rodents of the Caribbean: origin and diversification of hutias unravelled by nextgeneration museomics[END_REF]2017).

Disentangling the relationships of hutias is therefore crucial not only to define their highertaxonomic status but also to shed light on their potential biogeographical origin. Recent molecular phylogenies have suggested close affinities with the echimyid subfamily Euryzygomatomyinae (Clyomys, Euryzygomatomys, Trinomys;[START_REF] Fabre | Mitogenomic phylogeny, diversification, and biogeography of South American spiny rats[END_REF]) from Eastern Brazil. Molecular dating analyses have supported an Early to Mid-Miocene divergence of hutias [START_REF] Fabre | Rodents of the Caribbean: origin and diversification of hutias unravelled by nextgeneration museomics[END_REF]2017) with a suspected dispersal from an Amazonian ancestor [START_REF] Upham | Diversification and biogeography of the Neotropical caviomorph lineage Octodontoidea (Rodentia: Hystricognathi)[END_REF]. However, the recent sampling of a fast-evolving echimyid from Eastern Brazil (Carterodon sulcidens) has hindered our view of their basal relationships [START_REF] Upham | Evolution of caviomorph rodents: a complete phylogeny and timetree for living genera[END_REF][START_REF] Fabre | Mitogenomic phylogeny, diversification, and biogeography of South American spiny rats[END_REF].

After more than two decades spent on this complex clade, the higher-level phylogeny of Echimyidae still appears unresolved with a basal polytomy between Echimyidae subfamilies and Capromyidae, with several nodes poorly supported [START_REF] Fabre | Mitogenomic phylogeny, diversification, and biogeography of South American spiny rats[END_REF]). On the one hand, recent morphological studies managed to gather large datasets of cranio-dental, jaws and external characters [START_REF] Candela | Tooth morphology of Echimyidae (Rodentia, Caviomorpha): homology assessments, fossils, and evolution[END_REF][START_REF] Verzi | Phylogeny and Evolutionary Patterns of South American Octodontoid Rodents[END_REF]), but those were shown to be very susceptible to homoplasy [START_REF] Emmons | A Revision of the Genera of Arboreal Echimyidae (Rodentia: Echimyidae, Echimyinae), With Descriptions of Two New Genera[END_REF][START_REF] Candela | Tooth morphology of Echimyidae (Rodentia, Caviomorpha): homology assessments, fossils, and evolution[END_REF]. Molecular studies on the other hand mainly relied on few nuclear and mitochondrial genes [START_REF] Huchon | Rodent Phylogeny and a Timescale for the Evolution of Glires: Evidence from an Extensive Taxon Sampling Using Three Nuclear Genes[END_REF][START_REF] Galewski | Ecomorphological diversification among South American spiny rats (Rodentia; Echimyidae): a phylogenetic and chronological approach[END_REF][START_REF] Fabre | Diversification of South American spiny rats (Echimyidae): a multigene phylogenetic approach: A multigene molecular phylogeny of South American spiny rats[END_REF][START_REF] Stamatakis | RAxML version 8: a tool for phylogenetic analysis and post-analysis of large phylogenies[END_REF]2017;[START_REF] Álvarez | patterns and size evolution in caviomorph rodents[END_REF], along with a mitogenomic study [START_REF] Fabre | Mitogenomic phylogeny, diversification, and biogeography of South American spiny rats[END_REF]). However, the haploidy and maternal inheritance of mitochondrial DNA reduce its effective population size, leading to a fast fixation of mutations [START_REF] Brown | Rapid evolution of animal mitochondrial DNA[END_REF], which increases the odds of saturation and makes sequences less reliable to infer deep relationships [START_REF] Springer | Mitochondrial Versus Nuclear Gene Sequences in Deep-Level Mammalian Phylogeny Reconstruction[END_REF].

Moreover, the absence of recombination in the mitochondrial genome means that all mitochondrial genes belong to a single loci. This is problematic as it is known that gene trees may differ from the species tree [START_REF] Maddison | Gene Trees in Species Trees[END_REF][START_REF] Huang | Cause of gene tree discord? Distinguishing incomplete lineage sorting and lateral gene transfer in phylogenetics[END_REF]. These points advocate the use of several diverse nuclear DNA loci as new molecular markers to explore the Echimyidae phylogeny.

In this paper we screened the OrthoMaM public database [START_REF] Ranwez | OrthoMaM: A database of orthologous genomic markers for placental mammal phylogenetics[END_REF]Douzery et al. 2014) in order to gather a dataset of approximately 500 orthologous nuclear exons with contrasting molecular evolution properties in regard to their evolutionary rate, distribution of substitutions, and nucleotide content. This strategy would allow bypassing phylogenetic inference limitations induced by the analysis of a restricted class of marker (e.g. slow-evolving or G+C-rich markers). We then used an exon capture NGS protocol to amplify and sequence these coding genes markers for 35 non-model Echimyidae species, bringing together the largest molecular dataset to date for a rodent family. Bayesian and maximum-likelihood phylogenomics allowed us to generate the first fully-resolved higher-level Echimyidae phylogeny. We therefore explored three questions regarding Echimyidae relationships: (i) How do we resolve the deepest divergences among echimyid subfamilies and other previously unresolved nodes of their phylogeny? (ii) Can we settle the taxonomic status of the Capromyidae (hutias)? and (iii) How do the evolutionary rate, substitution distribution and nucleotide content of nuclear DNA exons impact phylogenetic inferences?

MATERIALS AND METHODS

Taxonomic sampling

Biological samples were obtained from both museum and fresh tissue collections from the Museum of Comparative Zoology (MCZ: Cambridge, MA), the Museum of Vertebrate Zoology, University of California (MVZ: Berkeley, CA), the U.S. National Museum of Natural History (USNMH-Smithsonian: Washington DC), the Nationaal Natuurhistorisch Museum (Naturalis:

Leiden, The Netherlands), the Royal Ontario Museum (ROM: Toronto, Canada), the Universidade Federal do Espírito Santo (UFES-CTA Animal Tissue Collection : Vitória, Brazil), and the University of Montpellier (UM-ISEM: Montpellier, FR), and are the same as those used in [START_REF] Fabre | Mitogenomic phylogeny, diversification, and biogeography of South American spiny rats[END_REF] for mitogenomics analyses. We re-identified the specimen of Mesomys stimulax RMNHN:MAM 21728 from Surinam, that was sequenced in [START_REF] Fabre | Mitogenomic phylogeny, diversification, and biogeography of South American spiny rats[END_REF]. This specime was identified by [START_REF] Husson | The mammals of Surinam[END_REF]p. 438-440) and sampled near the Vier Gebroeders mountain near Sipaliwini in the South Eastern part of Nickerie District. This specimen was misidentified in [START_REF] Fabre | Mitogenomic phylogeny, diversification, and biogeography of South American spiny rats[END_REF] as Mesomys stimulax, and we re-identified it as Mesomys hispidus. In total, we obtained tissues for 35 taxa spanning almost all genera of Echimyidae and Capromyidae. In addition, we searched NCBI for complete genomes of two Phiomorpha rodents (Heterocephalus glaber and Fukomys damarensis, Genbank assembly accession GCA_000247695.1 and GCA_000743615.1 respectively), and three Caviomorpha rodents Chinchilla lanigera (Chinchilloidea), Cavia porcellus (Cavioidea), and Octodon degus (Octodontidae), with accession numbers GCA_000276665.1, GCA_000151735.1 and GCA_000260255.1 respectively. These lineages span most Hystricognathi clades [START_REF] Upham | Evolution of caviomorph rodents: a complete phylogeny and timetree for living genera[END_REF] and served as outgroups in our phylogeny analyses. Due to lack of DNA or absence of tissues, five other known Echimyinae genera (Diplomys, Leiuromys, Lonchothrix, Santamartamys and Toromys) were not included in this study.

Exon sampling

Prior to the sequence capture, we identified three molecular evolution characteristics that may influence the quality and quantity of phylogenetic signal of each nuclear DNA exon: its evolutionary rate, its among-site substitution rate heterogeneity and its nucleotide composition [START_REF] Betancur-R R | Addressing gene tree discordance and non-stationarity to resolve a multi-locus phylogeny of the flatfishes (Teleostei: Pleuronectiformes)[END_REF][START_REF] Townsend | Profiling Phylogenetic Informativeness[END_REF][START_REF] Romiguier | Less Is More in Mammalian Phylogenomics: AT-Rich Genes Minimize Tree Conflicts and Unravel the Root of Placental Mammals[END_REF]. The aim of this study was to capture 500 exons combining various parameter values of these three variables that are known to impact phylogenetic analyses. To do so, we used the orthologous gene database OrthoMaM [START_REF] Douzery | OrthoMaM v8: A Database of Orthologous Exons and Coding Sequences for Comparative Genomics in Mammals[END_REF], which provides sequences and molecular evolution characteristics for mammalian nuclear exons. First, evolutionary rates of the exons at the hystricognath taxonomic scale were approximated by their Relative Rate of Evolution (RER) at the mammalian scale, as computed with a Super Distance Matrix approach and provided in OrthoMaM [START_REF] Ranwez | OrthoMaM: A database of orthologous genomic markers for placental mammal phylogenetics[END_REF]. Second, we used the shape parameter of the Gamma distribution that models the among-site heterogeneity in rates of evolution (α) as a proxy of variability distribution. When 0 < α < 1, the distribution is skewed towards a majority of invariable sites and a small number of positions accumulating most of the variability. In contrast when α > 1, substitutions are more evenly distributed among sites [START_REF] Yang | Among-site rate variation and its impact on phylogenetic analyses[END_REF]. Finally, nucleotide composition was estimated by the mean G+C content on the third codon positions of alignments (GC3) [START_REF] Clay | GC3 of Genes Can Be Used as a Proxy for Isochore Base Composition: A Reply to Elhaik et al[END_REF]. Two focal species were chosen as a starting point for this exon selection: Mus musculus is the most represented rodent in the OrthoMaM database, and Cavia porcellus is a caviomorph lineage more closely-related to the Echimyidae ingroup. Exons longer than 400 base pairs (bp) were classified in categories depending on these three molecular evolution parameters: four RER categories ranging from slower-to faster-evolving markers (i.e., RER < 1, 1≤ RER < 2, 2 ≤ RER < 3 and 3 ≤ RER < 5), four α categories ranging from strong to moderate among-site rate heterogeneity (i.e., 0.25 < α < 0.50, 0.50 ≤ α < 0.75, 0.75 ≤ α < 1, and 1 ≤ α ), as well as 61 GC3 categories ranging from 29% to 94% by steps of 1%. Among the 4 (RER) x 4 (α) x 61 (GC3) = 976 possible combinations, 483 exons matched the 3-parameter criteria and were finally identified and selected (see Table S1).

Exon capture and assembly

The hystricognath rodent Octodon degus is the closest relative of Echimyidae for which full genome data are available [START_REF] Upham | Evolution of caviomorph rodents: a complete phylogeny and timetree for living genera[END_REF]. We thus used this species' genome to design capture probes for our Echimyidae dataset. The sequences of the 483 target exons were extracted from O. degus and other outgroups genome with the TBLASTN software [START_REF] Camacho | BLAST+: architecture and applications[END_REF]) using amino-acid sequence of the Mus musculus and Cavia porcellus exons from OrthoMaM.

The parameters used (wordsize = 3 and evalue < e -10 ) allowed to recover slightly divergent sequences, and only the best hit was retained for analysis. We then designed a capture probe set matching the 483 O. degus exons with a tiling (probe coverage) of 3x, which comprised 19,017 120-nucleotide long baits. DNA libraries for 35 Echimyidae species were subsequently prepared according to [START_REF] Tilak | A costeffective straightforward protocol for shotgun Illumina libraries designed to assemble complete mitogenomes from non-model species[END_REF]. Exons were then captured using a MyBaits (MYcroarray, Ann Arbour, USA) kit and following manufacturer protocol. Captured DNA was sequenced on a Illumina HiSeq 2500 platform (MGX, Montpellier, France). Raw 125-nucleotide long sequencing reads were paired in the Geneious software [START_REF] Kearse | Geneious Basic: An integrated and extendable desktop software platform for the organization and analysis of sequence data[END_REF]. We used a Spearman rank correlation and an Analysis of Variance as implemented in the R software (R. [START_REF] Team | R: A language and environment for statistical computing[END_REF] to check if the age or location of tissues could explain the number of reads obtained. Adapters sequences were also trimmed and remaining reads were subsequently mapped onto O. degus reference sequences, while allowing for a maximum of 15% mismatch per read, and a consensus was produced. A second assembly run was performed using this consensus as reference while allowing only 5% mismatch in order to filter out spurious reads. After removal of duplicated reads, we excluded regions with coverage less than 3 from subsequent analyses.

Alignment and orthology assessment

Homologous DNA exon sequences were aligned with the MACSE software, which aligns coding DNA sequences using their amino-acid translation while allowing for potential frame shifts [START_REF] Ranwez | MACSE: Multiple Alignment of Coding SEquences Accounting for Frameshifts and Stop Codons[END_REF]. After this multiple alignment step, we devised the following procedure to remove the so-called "outlier" sequences, i.e., misaligned, paralogous or illegitimately-captured sequences. A maximum likelihood (ML) tree was inferred from each exon alignment using RAxML v8.2.8 [START_REF] Stamatakis | RAxML version 8: a tool for phylogenetic analysis and post-analysis of large phylogenies[END_REF]) under the GTR+Γ model of DNA sequence evolution [START_REF] Rodríguez | The general stochastic model of nucleotide substitution[END_REF][START_REF] Yang | Among-site rate variation and its impact on phylogenetic analyses[END_REF]. In each tree, patristic distances between each taxon and the ancestral node of the crown ingroup (all Echimyidae and Capromyidae) were computed. As "outlier" sequences are likely to result into a topological misplacement and longer branches leading to it, we thus analyzed variations in the patristic distances with R (R. Core Team 2018). We used a linear model that included two independent qualitative predictor variables: the exon (with 483 modalities) and the species (41 modalities) involved in the corresponding log-transformed root-to-tip distance. We then flagged as "outlier" every sequence involved in a root-to-tip distance with a residual value of the model greater than 3 standard deviations. Outlier sequences were then pruned out from corresponding alignments. After this curating step, sequences were re-aligned with MACSE using default parameters, and sites with more than 20% of ambiguous character states and gaps were removed using trimAl v1.3. [START_REF] Capella-Gutiérrez | trimAl: a tool for automated alignment trimming in large-scale phylogenetic analyses[END_REF].

The nucleotide and amino-acid datasets can be accessed in the following GitHub repository: https://github.com/CourcelleM/Courcelle_et_al2018_Echimyidae Zipped nexus files contain the raw and curated alignments of 483 exons for the 40 species studied.

Phylogenetic analyses

All curated exon alignments were concatenated leading to a supermatrix with 601,816 sites, comprising 14 % missing character states. At the protein level, the supermatrix comprised 194,088 amino acid sites, with 9% missing character states. Previous analyses highlighted short divergence times and fast molecular evolution of some Echimyidae clades [START_REF] Lara | The Simultaneous Diversification of South American Echimyid Rodents (Hystricognathi) Based on Complete Cytochrome b Sequences[END_REF][START_REF] Fabre | Diversification of South American spiny rats (Echimyidae): a multigene phylogenetic approach: A multigene molecular phylogeny of South American spiny rats[END_REF]2017). For example, the uncertain position of the long branch of Carterodon sulcidens destabilized the species tree of [START_REF] Fabre | Mitogenomic phylogeny, diversification, and biogeography of South American spiny rats[END_REF]. An evolutionary model able to efficiently deal with homoplasies and long-branch attraction was thus recommended in order to reconstruct the Echimyidae phylogeny. A Bayesian inference was conducted on the concatenated exon alignments, both at the nucleotide and amino-acid levels, using the CAT-GTR+Γ mixture model implemented in PhyloBayes-MPI v1.8 [START_REF] Lartillot | PhyloBayes 3: a Bayesian software package for phylogenetic reconstruction and molecular dating[END_REF]. This sequence evolution model takes across-site composition and substitution rates heterogeneity into account and is thus well suited to deal with long-branch attraction artifacts [START_REF] Lartillot | A Bayesian Mixture Model for Across-Site Heterogeneities in the Amino-Acid Replacement Process[END_REF]. Two Markov Chains Monte Carlo (MCMC) were run independently and parameters were sampled every 5 cycles. Chain convergence was assessed when the two chains inferred the same bipartitions frequencies. Node posterior probabilities were computed after respectively 1,210 and 1,245 cycles (corresponding to approximately 1,000,000 generations) with a 30% burn-in discarding the first 400 trees.

A potential limitation of the supermatrix approach is that the phylogenetic signal brought by each specific gene might be diluted during the joint analyses [START_REF] Song | Resolving conflict in eutherian mammal phylogeny using phylogenomics and the multispecies coalescent model[END_REF]. This can lead to incongruences when the genes bear conflicting information, for example due to incomplete lineage sorting or hidden paralogy [START_REF] Kubatko | Inconsistency of Phylogenetic Estimates from Concatenated Data under Coalescence[END_REF][START_REF] Kuraku | Impact of asymmetric gene repertoire between cyclostomes and gnathostomes[END_REF]. We thus also estimated the echimyid species tree by a multi-coalescent approach. To account for phylogenetic uncertainty, 100 bootstrap trees were inferred for each exon under the GTR+Γ model with RaxML. The resulting collection of 48,300 bootstrapped gene trees were then inputed to the ASTRAL-2 software [START_REF] Mirarab | ASTRAL-II: coalescent-based species tree estimation with many hundreds of taxa and thousands of genes[END_REF]. As "gene+site -wise" bootstrapping may generate artefactual conflicts among gene trees [START_REF] Simmons | Gene-wise resampling outperforms sitewise resampling in phylogenetic coalescence analyses[END_REF], we repeated this analysis using the "--gene-only" bootstrapping option of ASTRAL-III [START_REF] Zhang | ASTRAL-III: polynomial time species tree reconstruction from partially resolved gene trees[END_REF]). An assumption of summary coalescence methods like ASTRAL is that there is no intra-locus recombination. According to [START_REF] Scornavacca | Incomplete lineage sorting in mammalian phylogenetics[END_REF], the exon-level at which we analysed our dataset is likely the most robust for phylogenomic studies. To ascertain that intra-locus recombination did not distort our results, we re-analyzed our dataset using the SVDquartets method (Chifman and Kubatko 2015) as implemented in PAUP* [START_REF] Swofford | PAUP*: phylogenetic analysis using parsimony (* and other methods)[END_REF]. As SVDquartets uses quartets of taxa to infer phylogenetic relationships at the individual alignment site-level, it does not require the assumption that there is no intralocus recombination.

A corollary point relative to our analysis protocol is the impact level on phylogenetic inferences of the large spectrum of evolutionary rates, α values and nucleotide compositions of exons sampled (see [START_REF] Betancur-R R | Addressing gene tree discordance and non-stationarity to resolve a multi-locus phylogeny of the flatfishes (Teleostei: Pleuronectiformes)[END_REF][START_REF] Townsend | Profiling Phylogenetic Informativeness[END_REF], and Romiguier et al. 2013 respectively). We therefore measured the congruence between the reference supermatrix tree and phylogenetic trees inferred along a molecular evolution property gradient. Exons were sorted out according to their (1) rate of evolution, (2) α value and (3) nucleotide composition and were then concatenated. This ensured us to obtain three collections of aligned exons with ascending values of RER, α and GC3 respectively. We then used consecutive bins of 15,000 sites from each matrix to analyze datasets of the same length but with gradually contrasted molecular properties. We inferred 100 bootstrap trees for each bin and computed their distance to the "reference" tree using the Kscore metric, which takes into account the topological and branch length differences [START_REF] Soria-Carrasco | The K tree score: quantification of differences in the relative branch length and topology of phylogenetic trees[END_REF]). As we are interested in the K-score difference between the bins, the choice of the reference tree was not of prime importance, and we conventionally chose to use the phylogenetic tree yielded by our supermatrix Bayesian analysis.

Molecular dating

Molecular dating analyses were performed with PhyloBayes v4.1 [START_REF] Lartillot | PhyloBayes 3: a Bayesian software package for phylogenetic reconstruction and molecular dating[END_REF]). We used a Bayesian relaxed molecular clock with rate autocorrelation along the branches of the tree [START_REF] Thorne | Estimating the rate of evolution of the rate of molecular evolution[END_REF][START_REF] Thorne | Divergence Time and Evolutionary Rate Estimation with Multilocus Data[END_REF] on the reference topology yielded by the supermatrix analysis. Five uniform fossil constraints were employed to calibrate the node of the phylogeny, with soft bounds to allow dates to be drawn in an interval 5 % greater than the set boundaries [START_REF] Yang | Bayesian Estimation of Species Divergence Times Under a Molecular Clock Using Multiple Fossil Calibrations with Soft Bounds[END_REF]. The younger bound of the root was set to the oldest hystricognaths of South America at 41.0 Million years ago (Mya) [START_REF] Antoine | Middle Eocene rodents from Peruvian Amazonia reveal the pattern and timing of caviomorph origins and biogeography[END_REF]). The older bound was set to the early Eocene at 56.0 Mya. In the outgroup, the divergence time among Chinchilloidea and Octodontoidea was set between the oldest stem Chinchilloidea member †Eoviscaccia frassinettii (31.3 Mya, Bertrand et al. 2012) and the earliest caviomorph fossils (41.0 Mya, [START_REF] Antoine | Middle Eocene rodents from Peruvian Amazonia reveal the pattern and timing of caviomorph origins and biogeography[END_REF]). In the ingroup, the first calibration used was the divergence between Trinomys and Clyomys+Euryzygomatomys genera (Euryzygomatomyinae subfamily). Based on the †Theridomysops parvulus age, the oldest stem lineage to Euryzygomatomys+Clyomys genera, a lower bound was set to 6.0 Mya for the Euryzygomatomyinae node [START_REF] Verzi | Phylogeny and Evolutionary Patterns of South American Octodontoid Rodents[END_REF]). The upper bound was subsequently set to 11.8 Mya, the deepest level of the Laventan SALMA (South American Land Mammal Age) [START_REF] Van Der Geer | Phylogeny and adaptive diversity of rodents of the family Ctenomyidae (Caviomorpha): delimiting lineages and genera in the fossil record[END_REF]. We then constrained the most recent common ancestor of the Myocastorini to range between the oldest stem taxon to Thrichomys (6.0 Mya, [START_REF] Olivares | Phylogenetic affinities of the late Miocene echimyid †Pampamys and the age of Thrichomys (Rodentia, Hystricognathi)[END_REF][START_REF] Verzi | Phylogeny and Evolutionary Patterns of South American Octodontoid Rodents[END_REF]) and the start of the Laventan SALMA (11.8 Mya, [START_REF] Van Der Geer | Phylogeny and adaptive diversity of rodents of the family Ctenomyidae (Caviomorpha): delimiting lineages and genera in the fossil record[END_REF].

Finally, we dated the lower bound of the Echimyini divergence with the oldest representative of the Echimys+Phyllomys clade : †Maruchito trilofodonte (15.7 Mya, [START_REF] Vucetich | Los roedores de la Formación Collón Cura (Mioceno Medio) y la Ignimbrita Pilcaniyeu. Cañadón del Tordillo, Neuquén[END_REF][START_REF] Kay | Vertebrate Paleontology in the Neotropics: The Miocene Fauna of La Venta, Colombia[END_REF]). The upper bound was set to the start of the Deseadan SALMA (24.2 Mya, [START_REF] Dunn | A new chronology for middle Eocene-early Miocene South American Land Mammal Ages[END_REF].

Locomotor habits reconstruction

The evolution of locomotor habits of echimyids was modeled using the R package diversitree (FitzJohn 2012). Coding for locomotion lifestyle was as follows: (1) semi-arboreal, scansorial, or arboreal (Callistomys, Capromys, Dactylomys, Echimys, Isothrix, Kannabateomys, Mesocapromys, Makalata, Pattonomys, Plagiodontia, Mesomys, Mysateles, Olallamys, Phyllomys); (2) semifossorial (Clyomys, Euryzygomatomys, Carterodon); (3) terrestrial (Proechimys, Geocapromys, Hoplomys, Thrichomys, Trinomys); and (4) semi-aquatic (Myocastor). Ancestral character state estimations were implemented in the diversitree package in the R programming language (R. Core Team 2018) using the "make.mkn" function. This maximum likelihood (ML) mapping approach fits the model of [START_REF] Pagel | Detecting correlated evolution on phylogenies: a general method for the comparative analysis of discrete characters[END_REF] to reconstruct states at all ancestral nodes of our phylogeny using the Mk2/Mk-n models. Three models of trait evolution were fitted and compared: (i) all rates for locomotor habit transitions are equal, (ii) rates are symmetrical, and (iii) all rates are allowed to be different. The first model provided the lowest AIC. Plotting of the results were performed using a R function from [START_REF] Moen | Testing convergence Versus History: Convergence Dominates Phenotypic Evolution for over 150 Million Years in Frogs[END_REF].

RESULTS AND DISCUSSION

Exon capture as an efficient tool for phylogenomic analysis

In this study, we performed an exon capture experiment in order to collect a vast molecular dataset for species for which previous analyses were limited to mitochondrial genomes and a handful of nuclear exons [START_REF] Lara | The Simultaneous Diversification of South American Echimyid Rodents (Hystricognathi) Based on Complete Cytochrome b Sequences[END_REF][START_REF] Huchon | From the Old World to the New World: A Molecular Chronicle of the Phylogeny and Biogeography of Hystricognath Rodents[END_REF][START_REF] Leite | Evolution of South American spiny rats (Rodentia, Echimyidae): the star-phylogeny hypothesis revisited[END_REF][START_REF] Upham | Evolution of caviomorph rodents: a complete phylogeny and timetree for living genera[END_REF][START_REF] Fabre | Mitogenomic phylogeny, diversification, and biogeography of South American spiny rats[END_REF][START_REF] Álvarez | patterns and size evolution in caviomorph rodents[END_REF]. Given a reference sequence, exon capture allows to easily isolate and sequence DNA from other individuals from the same species, but its prime interest is to extend this capacity to genes from sister-group taxa that are likely to exhibit more sequence divergence. If a close enough reference is available for non-model species, this technique is a viable alternative to shotgun sequencing for many investigations where molecular data is required, from ancient DNA studies to whole-genome sequencing and phylogenomics [START_REF] Enk | Ancient Whole Genome Enrichment Using Baits Built from Modern DNA[END_REF][START_REF] Mccormack | Sequence capture of ultraconserved elements from bird museum specimens[END_REF][START_REF] Teasdale | Identification and qualification of 500 nuclear, single-copy, orthologous genes for the Eupulmonata (Gastropoda) using transcriptome sequencing and exon capture[END_REF]. The capture approach is limited by the raw genetic distance between reference and target sequences: hybridization efficiency decreases when probes increasingly differ from target sequences, resulting into reduced efficiency of captures and amplifications. The protocol followed in the present study allowed us to capture exon sequences from Echimyidae and Capromyidae using an octodontid representative (here Octodon degus), for sequence divergence up to 20% (raw DNA pairwise distance from O. degus to other species, computed from final exon alignments). This is in agreement with other studies that reported capture enrichment success for divergent sequences at least up to 15% [START_REF] Bragg | Exon capture phylogenomics: efficacy across scales of divergence[END_REF], and up to 39% with modified hybridization protocols [START_REF] Li | Capturing protein-coding genes across highly divergent species[END_REF]. We thus validate this exon capture protocol for two families whose divergence time is estimated around 23.6 Mya [START_REF] Upham | Evolution of caviomorph rodents: a complete phylogeny and timetree for living genera[END_REF], 31 Mya (Álvarez et al. 2017), or 21.1

Mya ([25.1 -18.1] Mya in the present analysis). The number of paired raw reads captured for each of the 35 echimyid taxa sequenced here ranged from 119,544 to 10.4 Million (Table 1). Interspecies read number heterogeneity is often explained by a lower quality of the biological material because degraded DNA is more difficult to sequence and assemble through new generation sequencing [START_REF] Burrell | The use of museum specimens with high-throughput DNA sequencers[END_REF]. This is not the case here, as the number of raw reads was not correlated to the age of tissues (Spearman's rank correlation coefficient, p-value = 0.4), nor explained by the origin of the sample tissue (ANalysis Of VAriance, p-value = 0.16), two factors known to hamper NGS results [START_REF] Almeida | Each flying fox on its own branch: A phylogenetic tree for Pteropus and related genera (Chiroptera: Pteropodidae)[END_REF]. Despite variation in the number of reads available for each taxon, 89 to 99% of the 483 targeted exons were captured for each species (Table 1). Moreover, our methodology allowed us to capture a sample of nuclear exons covering a wide range of molecular and evolutionary characteristics (Fig. 1). The relative evolutionary rates (RER) of exons computed from the gene trees inferences ranged between approximately 0.01 (very conserved) and 6 (fast-evolving), denoting a 600-fold difference in the exons rates of evolution. The α parameter varies from 0.02 to 13.01. The highest value corresponds to the first exon of the gene TEX10 (Ensembl reference ENSG00000136891). It was a clear outlier as the second highest α was estimated at approximately 2.2. Checking the alignment of TEX10 revealed that the sequence is extremely conserved (RER of 0.01), but rare substitutions seemed to occur at random third codon positions along the sequence, leading to an almost uniform probability of substitution for each site (very high α value). It is interesting to note that the distribution of α and GC3 values show respectively 82% and 86% overlap between mammalian and rodent scales (Fig. 1B and1C). It thus seems that the estimates of α and GC3 values do not significantly depend upon the taxonomic scale at which they are computed. Values from the OrthoMam database therefore seem to be a reasonable approximation for mammalian species at a lower taxonomic scale. Relative evolution rates however are found generally lower in our echimyids dataset than in OrthoMam (63% distribution overlap, Fig. 1A). Since our dataset contains phylogenetically closer taxa than those in the Orthomam database, exons show less variability among the species studied here, which leads to lower relative evolutionary rates.

Even if our results highlight a broad range of exon capture success, we did observe a negative correlation between assembly coverage and the relative rate of evolution of captured exons (Fig. 2B). This may be a hint that divergent sequences are harder to capture: the more they differ from the O. degus reference, the less efficient is the hybridization between probes and target sequences (Fig. 2A and2B). Moreover, optimal exon coverage was expected for a balanced nucleotide content, i.e. a GC3 content of 50% (see e.g. [START_REF] Bi | Transcriptome-based exon capture enables highly cost-effective comparative genomic data collection at moderate evolutionary scales[END_REF]. Indeed, a strong nucleotide bias lowers sequence complexity, and therefore the probes specificity. Additionally, amplification steps are less efficient when sequences are richer in G and C nucleotides [START_REF] Viswanathan | Template secondary structure promotes polymerase jumping during PCR amplification[END_REF].

In this study, mean exon coverage strictly decreased with GC3 content (Fig. 2C). This relation suggests that the lower complexity of probes in AT-rich regions may be less detrimental than the difficulties brought by G+C-rich sequences, such as stronger DNA secondary structure hampering amplification steps [START_REF] Kozarewa | Amplification-free Illumina sequencing-library preparation facilitates improved mapping and assembly of (G+C)-biased genomes[END_REF], or higher probability of sequencing error [START_REF] Nakamura | Sequence-specific error profile of Illumina sequencers[END_REF]. Thus, the A+T-richness of sequences does not seem to have been a limiting factor to capture exon at the rodent family taxonomic scale, as illustrated by the higher sequencing depth of A+T-rich exons (Fig. 2C).

A fully resolved Echimyidae phylogeny

The phylogeny of South American spiny rats greatly benefited from the advent of highthroughput sequencing and molecular phylogeny. Morphological analyses [START_REF] Carvalho | Relationships among extant and fossil echimyids (Rodentia: Hystricognathi)[END_REF][START_REF] Emmons | A Revision of the Genera of Arboreal Echimyidae (Rodentia: Echimyidae, Echimyinae), With Descriptions of Two New Genera[END_REF]) and molecular star-phylogenies [START_REF] Lara | The Simultaneous Diversification of South American Echimyid Rodents (Hystricognathi) Based on Complete Cytochrome b Sequences[END_REF][START_REF] Leite | Evolution of South American spiny rats (Rodentia, Echimyidae): the star-phylogeny hypothesis revisited[END_REF] were followed by probabilistic methods using more sophisticated evolutionary models which, even if applied to only a handful of mitochondrial and nuclear genes, allowed a partial resolution of the Echimyidae phylogeny [START_REF] Galewski | Ecomorphological diversification among South American spiny rats (Rodentia; Echimyidae): a phylogenetic and chronological approach[END_REF][START_REF] Upham | Evolution of caviomorph rodents: a complete phylogeny and timetree for living genera[END_REF][START_REF] Fabre | Diversification of South American spiny rats (Echimyidae): a multigene phylogenetic approach: A multigene molecular phylogeny of South American spiny rats[END_REF]2017;[START_REF] Álvarez | patterns and size evolution in caviomorph rodents[END_REF][START_REF] Emmons | A Review of the Pattonomys/Toromys Clade (Rodentia: Echimyidae), with Descriptions of a New Toromys Species and a New Genus[END_REF]. Still, several nodes were poorly supported or did not gain maximal statistical support [START_REF] Fabre | Mitogenomic phylogeny, diversification, and biogeography of South American spiny rats[END_REF], despite the use of complete mitogenomes and up to five nuclear exons. Here, two independent phylogenetic approaches were used on our large nuclear dataset to untangle the Echimyidae higher-level phylogeny: (1) a Bayesian inference on the supermatrix of all exon sequences and (2) multicoalescent-based approaches using the gene trees as a primary input. The two analyses yielded highly similar results (Fig. 3): the same topology, and each node unambiguously supported, with respectively posterior probabilities and local posterior probabilities values equals to 1 (Fig. 3). The tree obtained is the same regardless of the bootstrapping approach used ("gene+site -wise" or "gene-wise", see section 2.5). The SVDquartets analysis also produced the same topology. Thus, we here report the first fully and well-supported genus-level Echimyidae phylogeny.

The paraphyly of Echimyidae as currently understood is evidenced as it includes the Capromyidae as sister group to the Brazilian echimyid Carterodon. We also recovered three major clades : (1) Capromyinae (new rank), comprising living hutias from the West Indies (including Plagiodontia, Geocapromys, Capromys, Mysateles and Mesocapromys), (2) Euryzygomatomyinae including terrestrial and semifossorial taxa from eastern and central South America (Trinomys, Clyomys and Euryzygomatomys), and (3) Echimyinae, comprising the remaining, mostly arboreal spiny rats (Dactylomys, Echimys, Isothrix, Kannabateomys, Makalata, Mesomys, Ollalamys, Pattonomys, Phyllomys, Callistomys, Hoplomys, Myocastor, Proechimys and Thrichomys). These findings reinforce phylogenetic hypotheses proposed in previous studies [START_REF] Galewski | Ecomorphological diversification among South American spiny rats (Rodentia; Echimyidae): a phylogenetic and chronological approach[END_REF][START_REF] Upham | Evolution of caviomorph rodents: a complete phylogeny and timetree for living genera[END_REF][START_REF] Fabre | Mitogenomic phylogeny, diversification, and biogeography of South American spiny rats[END_REF], however the relative positions of the three subfamilies were conflictual. In line with most recent phylogenies [START_REF] Fabre | Rodents of the Caribbean: origin and diversification of hutias unravelled by nextgeneration museomics[END_REF]2017), our results firmly corroborate sister-relationships between Capromyinae and Carterodon (for which we create a new subfamily, see below), forming a clade sister to Euryzygomatomyinae. This clade is in turn sister to Echimyinae (Fig. 3).

Within Echimyinae, we again retrieved the surprising affinity between the semi-aquatic

Myocastor coypus and the arboreal Callistomys pictus [START_REF] Loss | Unexpected phylogenetic relationships of the painted tree rat Callistomys pictus (Rodentia: Echimyidae)[END_REF][START_REF] Fabre | Mitogenomic phylogeny, diversification, and biogeography of South American spiny rats[END_REF]. These taxa form a well-supported clade together with Thrichomys and Hoplomys + Proechimys, named Myocastorini sensu [START_REF] Fabre | Mitogenomic phylogeny, diversification, and biogeography of South American spiny rats[END_REF]. Other arboreal genera of the subfamily are grouped in the Echimyini tribe. The lineage leading to the brush-tailed rat Isothrix diverged first, followed by the lineage leading to Mesomys. This finding is at odds with most previous studies that proposed a clade Isothrix + Mesomys [START_REF] Fabre | Mitogenomic phylogeny, diversification, and biogeography of South American spiny rats[END_REF][START_REF] Álvarez | patterns and size evolution in caviomorph rodents[END_REF]. We recovered a rather short branch between these two genera (approximately 0.002 substitutions per site) that could not be reliably inferred using only the small gene sampling of previous studies [START_REF] Upham | Evolution of caviomorph rodents: a complete phylogeny and timetree for living genera[END_REF][START_REF] Fabre | Mitogenomic phylogeny, diversification, and biogeography of South American spiny rats[END_REF]. The bamboo rat clade is also retrieved with Kannabateomys amblyonyx being a sister-taxon to a clade formed by Dactylomys dactylinus and Olallamys albicauda as in previous studies of [START_REF] Galewski | Ecomorphological diversification among South American spiny rats (Rodentia; Echimyidae): a phylogenetic and chronological approach[END_REF] and [START_REF] Upham | Diversification and biogeography of the Neotropical caviomorph lineage Octodontoidea (Rodentia: Hystricognathi)[END_REF]. Species of Echimys are closely related to Phyllomys, and this clade is sister to Makalata, forming another clade sister to Pattonomys in accordance with recent phylogenetic results [START_REF] Fabre | Mitogenomic phylogeny, diversification, and biogeography of South American spiny rats[END_REF][START_REF] Emmons | A Review of the Pattonomys/Toromys Clade (Rodentia: Echimyidae), with Descriptions of a New Toromys Species and a New Genus[END_REF]. The subfamily Euryzygomatomyinae comprises Trinomys and its semifossorial sister-clade Clyomys laticeps and Euryzygomatomys spinosus. The Capromyinae includes Capromys, Mesocapromys, and Mysateles, which form the Capromyini tribe. This clade is the sister group to the Hispaniolan hutia Plagiodontia aedium (Plagiodontini).

As our phylogenetic analyses were performed on a wide sampling of nuclear exons with contrasting molecular properties, we evaluated whether the evolutionary rate, the substitution distribution and the G+C-content of exons distorted gene tree inferences. The variation of the congruence of trees reconstructed along a gradient of these three properties has been measured by the K-score metric (Fig. 4). We did not detect any clear impact of the molecular evolutionary rate, the value of the α parameter or the GC content on the level of the individual exon tree congruence.

We however observed an edge effect where trees inferred from the subsamples with (1) low Relative Evolution Rates (RER) and (2) low α values show higher K-score values than the average (Fig. 4). This suggests that the slowest-evolving exons, as well as those where multiple substitutions occur, are more likely to generate homoplasy. A closer examination of the phylogenetic trees revealed that rather than providing strongly incongruent relationships among Echimyidae, these samples lack resolution to resolve the shortest branches, such as the deep branching of the different subfamilies or the branches leading to Carterodon sulcidens or Mesomys / Lonchothrix (Fig. 3). As these issues are compensated by other exons in our dataset, we postulate that including contrasting molecular characteristics did not mislead our phylogenetic inferences but instead allowed us to reach such a high statistical resolution.

Taxonomic implications for the Capromyidae and Carterodon sulcidens

With the resolution of the phylogenetic relationships among the Echimyidae clades, this study brings the first unequivocal material to disentangle the higher-level echimyid classification.

First, the relationships of the Caribbean hutias, the Capromyidae, in regard to Echimyidae have been uncertain for a long time (Wood and Howland 1979). Molecular studies yielded two alternative topologies that question the taxonomic status of Capromyidae as well as the monophyly of Echimyidae (see [START_REF] Fabre | Mitogenomic phylogeny, diversification, and biogeography of South American spiny rats[END_REF]). The first hypothesis was the reciprocal monophyly of Capromyidae and Echimyidae [START_REF] Upham | Diversification and biogeography of the Neotropical caviomorph lineage Octodontoidea (Rodentia: Hystricognathi)[END_REF]. This view strengthened the status of Capromyidae as a distinct family, closest relative of extant echimyids. The second hypothesis recovered in the literature consisted in grouping Capromyidae and Euryzygomatomyinae within Echimyidae [START_REF] Galewski | Ecomorphological diversification among South American spiny rats (Rodentia; Echimyidae): a phylogenetic and chronological approach[END_REF][START_REF] Fabre | Diversification of South American spiny rats (Echimyidae): a multigene phylogenetic approach: A multigene molecular phylogeny of South American spiny rats[END_REF][START_REF] Stamatakis | RAxML version 8: a tool for phylogenetic analysis and post-analysis of large phylogenies[END_REF]2016). Echimyidae would here become paraphyletic, unless Capromyidae are actually viewed as a subfamily within Echimyidae.

Our results provide unambiguous support to the second hypothesis, i.e. a Capromyinae + Euryzygomatomyinae clade nested within Echimyidae (Fig. 3, Fig. 5). In the present study, we used both a wide taxonomic sampling, including all extant genera of Capromyidae, as well as the best character sampling of Echimyidae and Capromyidae rodents to date with more than 600,000 nuclear DNA and 200,000 protein sites. We also used slow-evolving exons with various molecular evolution characteristics (see previous parts) and different phylogenetic methodologies to decrease the impact of well-known methodological artifacts such as site saturation, long-branch artifacts [START_REF] Lartillot | A Bayesian Mixture Model for Across-Site Heterogeneities in the Amino-Acid Replacement Process[END_REF][START_REF] Delsuc | Phylogenomics and the reconstruction of the tree of life[END_REF], as well as incomplete lineage sorting and hidden paralogy [START_REF] Kubatko | Inconsistency of Phylogenetic Estimates from Concatenated Data under Coalescence[END_REF][START_REF] Kuraku | Impact of asymmetric gene repertoire between cyclostomes and gnathostomes[END_REF]. To fit our results as well as several recent phylogenetic inferences [START_REF] Fabre | Rodents of the Caribbean: origin and diversification of hutias unravelled by nextgeneration museomics[END_REF]2017) in a systematic framework, we thus propose to rank Capromyidae, as defined by [START_REF] Smith | Introduction to Mammals[END_REF], as a subfamily -Capromyinae -within Echimyidae definition sensu [START_REF] Mckenna | Classification of mammals: above the species level[END_REF]. Moreover, previous studies struggled to assert the phylogenetic placement of the enigmatic Carterodon sulcidens (but see also [START_REF] Vucetich | Un nuevo Echimyidae (Rodentia, Hystricognathi) de la Edad Colhuehuapense de Patagonia y consideraciones sobre la sistemática de la familia[END_REF], also referred to as incertae sedis [START_REF] Fabre | Family Echimyidae (Hutias, Coypu, and South American Spiny rats)[END_REF]2017). We actually provide here the first highly supported relationship of this semifossorial Brazilian species within the echimyid phylogeny. This finding is inconsistent with previous studies where C. sulcidens was more closely related to Euryzygomatomyinae [START_REF] Emmons | A Review of the Pattonomys/Toromys Clade (Rodentia: Echimyidae), with Descriptions of a New Toromys Species and a New Genus[END_REF][START_REF] Upham | Evolution of caviomorph rodents: a complete phylogeny and timetree for living genera[END_REF][START_REF] Álvarez | patterns and size evolution in caviomorph rodents[END_REF]).

Here C. sulcidens displays the longest isolated branch from the phylogeny, and its path to the Echimyidae ancestral node is 20 % longer than the mean for other species (Fig. 3). Its new relationship implies major biogeographical consequences for the Greater Antilleans hutias as well as convergent evolution of fossoriality in echimyids (see also [START_REF] Carvalho | Relationships among extant and fossil echimyids (Rodentia: Hystricognathi)[END_REF][START_REF] Candela | Tooth morphology of Echimyidae (Rodentia, Caviomorpha): homology assessments, fossils, and evolution[END_REF][START_REF] Fabre | Diversification of South American spiny rats (Echimyidae): a multigene phylogenetic approach: A multigene molecular phylogeny of South American spiny rats[END_REF]). Indeed, it has been shown that ecological shifts were common in Echimyidae, leading to unexpected relationships such as the clade formed by the arboreal Callistomys and the semi-aquatic Myocastor in a mostly terrestrial tribe [START_REF] Loss | Unexpected phylogenetic relationships of the painted tree rat Callistomys pictus (Rodentia: Echimyidae)[END_REF][START_REF] Fabre | Mitogenomic phylogeny, diversification, and biogeography of South American spiny rats[END_REF]. The phylogenetic position of Carterodon sulcidens implies that the semifossorial adaptations likely evolved at least twice in Echimyidae during the Middle-Miocene (Table 2, Figure S2): once along the branch leading to the Euryzygomatomys + Clyomys clade, and once along the branch of Carterodon. This lifestyle has been suggested to be an adaptation in terms of thermoregulation as well as a way to escape predators in the open environments of Cerrado and Caatinga (de Albuquerque et al. 2012), the South-American biomes which host all semifossorial Echimyidae species and several other semifossorial caviomorph rodents [START_REF] Fabre | Mitogenomic phylogeny, diversification, and biogeography of South American spiny rats[END_REF]. This lifestyle convergence with its related morphological adaptations may have misled early taxonomists in grouping together Carterodon sulcidens and semifossorial Euryzygomatomyinae on the basis of morphological characters [START_REF] Carvalho | Relationships among extant and fossil echimyids (Rodentia: Hystricognathi)[END_REF][START_REF] Emmons | A Revision of the Genera of Arboreal Echimyidae (Rodentia: Echimyidae, Echimyinae), With Descriptions of Two New Genera[END_REF][START_REF] Candela | Tooth morphology of Echimyidae (Rodentia, Caviomorpha): homology assessments, fossils, and evolution[END_REF].

However, some morphologists have already pointed the distinctiveness of Carterodon with its unique grooved upper incisors (Fig. 3, insert A) [START_REF] Bishop | An annotated list of Caviomorph rodents collected in North-Eastern Mato Grosso, Brazil[END_REF], and some autapomorphic cranio-mandibular characters [START_REF] Vucetich | Un nuevo Echimyidae (Rodentia, Hystricognathi) de la Edad Colhuehuapense de Patagonia y consideraciones sobre la sistemática de la familia[END_REF]. Taking together these molecular, biogeographical and morphological evidences, we propose to define a new subfamily Carterodontinae, which would comprise the extant monotypic genus Carterodon (Fig. 3 2015); canal of the infraorbital nerve in the infraorbital foramen absent [START_REF] Carvalho | Relationships among extant and fossil echimyids (Rodentia: Hystricognathi)[END_REF].

Distribution: central and western part of Brazil in the Cerrado biome, ranging from the states of Minas Gerais to Mato Grosso [START_REF] Bezerra | A review of the distribution, morphometrics, and habit of Owl's Spiny Rat Carterodon sulcidens (Lund, 1841)(Rodentia: Echimyidae)[END_REF].

Based on the results of this study, we propose a revised classification of living echimyid rodents above the species level.

3.4.Biogeographical and systematic perspectives

The Echimyidae phylogeny inferred from ca. 500 exons also bolsters the current hypothesis of the echimyid dispersal, with the Caribbean Capromyinae diverging from Neotropical lineages during the mid-Miocene [START_REF] Fabre | Rodents of the Caribbean: origin and diversification of hutias unravelled by nextgeneration museomics[END_REF]2016). Mid-Miocene is marked by the opening of the habitats in South America which have shaped Echimyidae diversification along with the emergence of open habitat specialists such as Carterodon, Clyomys, and Thrichomys [START_REF] Galewski | Ecomorphological diversification among South American spiny rats (Rodentia; Echimyidae): a phylogenetic and chronological approach[END_REF][START_REF] Fabre | Diversification of South American spiny rats (Echimyidae): a multigene phylogenetic approach: A multigene molecular phylogeny of South American spiny rats[END_REF]2016). This period of echimyid ecological and morphological diversifications is also accompanied by the diversification of the hutias, thanks to their arrival in the Greater Antilles.

This colonization seems to have occurred in the mid-Miocene as inferred by our molecular dating analysis (Fig. 5, Table 2) as well as prior dating of this subfamily [START_REF] Fabre | Diversification of South American spiny rats (Echimyidae): a multigene phylogenetic approach: A multigene molecular phylogeny of South American spiny rats[END_REF][START_REF] Stamatakis | RAxML version 8: a tool for phylogenetic analysis and post-analysis of large phylogenies[END_REF][START_REF] Upham | Evolution of caviomorph rodents: a complete phylogeny and timetree for living genera[END_REF]. The close affinity of hutias with Carterodon might imply a potential origin in central South America. Unnamed capromyids have been reported from Quebrada Honda, a middle Miocene site in Bolivia [START_REF] Mckenna | Classification of mammals: above the species level[END_REF], but none of the rodents from this site has been described in detail, and recent changes in taxonomy indicate that capromyids do not appear to be present there [START_REF] Croft | The middle Miocene (Laventan) Quebrada Honda fauna, southern Bolivia and a description of its notoungulates[END_REF]. Extinction has probably hindered the original continental source pool of hutia ancestors, letting Carterodon as the sole echimyid survivor of their ancient radiation. An interesting avenue is that several echimyid subfossils (e.g. Heteropsomyinae) are known from the Greater Antilles deposits ( [START_REF] Van Der Geer | Phylogeny and adaptive diversity of rodents of the family Ctenomyidae (Caviomorpha): delimiting lineages and genera in the fossil record[END_REF]. Some of these lineages (Boromys, Brotomys) share with Carterodon a simple occlusal molar pattern [START_REF] Nowak | Walker's Mammals of the World[END_REF]. Further investigations of morphological and molecular systematic of Carterodon as well as Caribbean rodent subfossils will be paramount. Such analysis will help us estimate the number of colonizations to the West Indies, as well as shed light on the extinct morphological disparity of the diverse Echimyidae. Fig. 4 K-score metric between the reference phylogenetic tree of Echimyidae (Fig. 3) and 40 trees inferred from 15,000 sites-long subsamples of our dataset with a gradient of Relative Rate of Evolution (A), a gradient of among-site rate variation α (B) and a gradient of G+C content on third codon positions (C). X-axes for these figures are not linear scales. Bins boxes are arbitrarily evenly spaced. Parameter values for each five bins are annotated. Values range approximately from 0.36 to 2.6 for the RER, 0.064 to 0.992 for α and 0.34 to 0.67 for the GC3.

  Both the supermatrix and multi-coalescent analyses place C. sulcidens as a deep-branching sistergroup to the Capromyinae. The clade formed by the Carterodontinae and the Capromyinae is defined by 175 exclusive amino-acid synapomorphies. Examples of exclusive synapomorphies with no missing taxa in our dataset include i) a Valine to Glutamic acid substitution at position homologous to site 887 of the human Protocadherin Fat 1 protein (UniprotKB: Q14517); ii) a Leucine to Serine substitution at position homologous to site 3009 of the human Apolipoprotein B protein (UniprotKB: P04114).

  Carterodon Waterhouse, 1848. Included genera: only the monotypic Carterodon, represented by Carterodon sulcidens (Lund, 1841). Diagnosis: small (<200g adult body weight, Bezerra et al. 2011), semifossorial rodent with short tail (ca. 50% head and body length); anterior face of the upper incisor grooved, protoloph absent in upper molariforms (Patton and Emmons

Fig. 2

 2 Fig. 2 Variation in sequencing coverage (reads per site) of the 483 captured nuclear exons. Median coverage was computed for each exon from the sequencing depth of the 36 Echimyidae species. (A) Coverage of exons as a function of the mean distance from O. degus to Echimyidae species in the corresponding alignment. (B) Coverage of exons according to the exon's relative rate of evolution. (C) Coverage of exons according to the mean G+C content in the 3rd codon position of the corresponding alignment.

Fig. 3

 3 Fig. 3 Phylogenetic tree of Echimyidae species inferred from a Bayesian analysis of the concatenation of 483 nuclear exons (601,816 DNA sites). Sequence evolution was described under the CAT + GTR + Γ mixture model. Posterior probabilities are reported for nodes of the tree. The phylogeny is rooted with the Octodontidae Octodon degus. For readability, other Hystricognathi outgroups have been pruned. Branch lengths are expressed in expected number of substitutions per site. Colors delimiting clades are given along the topology. (a) represents a ventral view of the skull of Carterodon sulcidens. The arrow highlights the distinctive grooved incisors of this species.

Fig. 5

 5 Fig. 5 Molecular dating of the Echimyidae family divergence. The topology of the chronogram was constrained to the result of our Bayesian analysis of the concatenation of 483 nuclear exons (Fig. 3). Green-colored bars represent the 95% credibility time interval inferred for nodes. The gray overlay displays the expected time window of the West Indies colonization by the Capromyinae subfamily. It has been set between the upper bound of node C (origin of Carterodontinae) and lower bound of the node D (diversification of Capromyinae). The timescale is expressed in million years.

Table 1

 1 Capture and sequencing statistics of the 483 nuclear exons for the 35 studied Echimyidae species

	Species	Raw reads	Unique reads	Unique reads percentage	Captured exons
	Callistomys pictus	10 446 228	1 258 230	12.04	477
	Capromys pilorides	1 302 422	269 201	20.67	461
	Carterodon sulcidens	6 165 400	1 466 951	23.79	476
	Clyomys laticeps	6 828 172	1 657 105	24.27	473
	Dactylomys dactylinus	3 233 010	1 118 435	34.59	462
	Echimys chrysurus	119 544	33 143	27.72	430
	Echimys saturnus	5 545 284	984 119	17.75	468
	Euryzygomatomys spinosus	7 079 304	1 723 493	24.35	433
	Geocapromys thoracatus	7 927 390	2 013 689	25.40	477
	Hoplomys gymnurus	8 108 526	1 926 392	23.76	472
	Isothrix pagurus	5 915 678	1 958 259	33.27	474
	Isothrix sinnamariensis	9 330 248	1 595 364	17.10	473
	Kannabateomys amblyonyx	6 632 906	1 293 378	19.50	472
	Makalata didelphoides	3 587 512	910 757	25.39	472
	Mesocapromys melanurus	3 220 430	768 595	23.87	478
	Mesocapromys nana	7 632 848	1 972 264	25.84	476
	Mesomys hispidus	9 871 930	1 983 491	20.09	478
	Myocastor coypus	364 238	142 555	39.14	459
	Mysateles prehensilis	1 536 955	646 248	42.05	475
	Olallamys albicauda	10 053 242	1 923 423	19.13	468
	Pattonomys carrikeri	7 661 732	1 732 204	22.61	475
	Phyllomys blainvillii	3 233 902	583 265	18.04	467
	Phyllomys dasythrix	5 140 752	991 961	19.30	475
	Phyllomys lundi	8 579 706	1 941 857	22.63	479
	Phyllomys mantiqueirensis	4 427 564	802 854	18.13	478
	Phyllomys pattoni	3 201 342	654 229	20.44	476
	Plagiodontia aedium	690 040	555 914	80.56	474
	Proechimys cuvieri	7 405 714	1 875 484	25.32	474
	Thrichomys apereoides	9 789 590	1 793 456	18.32	437
	Trinomys albispinus	168 264	33 830	20.11

Table 2

 2 Molecular dating of deep Echimyidae phylogenetic nodes in this analysis as compared to recent studies. Node labels are reported on Fig.5. Node ages and 95% credibility intervals are given in million years ago. Variation of molecular properties and evolution characteristics of the 483 target nuclear exons, as captured in this study (histogram) at the hystricognath rodents scale and according to the Orthomam reference at the mammals scale (dashed line distribution). (A) Relative evolutionary rates of exons, as computed from the RaxML trees using a Super Distance Matrix approach. Values range from 0.01 to 6.4. (B) Value of the shape parameter α of the Gamma distribution for the among-site substitution rate heterogeneity, estimated by maximum-likelihood during the inference of RaxML trees for each exon alignment. Values range from 0.02 to 2.2. One outlier data point corresponding to the first exon of the gene TEX10 (Ensembl reference ENSG00000136891) is not reported (see Results and Discussion). (C) Mean G+C content in the 3rd codon position of exon's alignments. Values range from 0.24 to 0.96.

	Node	This study	Fabre et al. 2017	Álvarez et al. 2017 Upham and Patterson 2015
	A	14.4 [18.8-13.2] 13.6	26.73	18.2
	B	14.1 [18.2-12.9] 12.7	25.36	17.4
	C	13.4 [16.1-12.3] NA	NA	NA
	D	8.6 [10.1-7.5]	7.6	15.56	9.8
	E	12.8 [14.0-11.7] 12.7	23.65	15.8
	F	13.1 [17.1-12.0] 13.3	25.05	17.1
	G	11.8 [15.7-10.7] 12.0	23.26	16.2
	H	10.8 [11.9-9.8] 9.4	22.49	15.5

Acknowledgments

We thank Fabien Condamine, François Catzeflis, Frédéric Delsuc, Nicolas Galtier, Vincent Ranwez, Nathan Upham and two anonymous referees for discussions about analyses and the redaction of this article. We are grateful to the following people and institutions for granting access to study skins and samples: Darrin Lunde, Nicole Edmison, and Louise Emmons (NMH, Smithsonian Institution), Chris Conroy and Jim Patton (MVZ Berkeley), Ana Carolina Loss and Leonora Costa (UFES), Judy Chupasko and Hopi Hoekstra at the Museum of Comparative Zoology (MCZ, Cambridge USA), Steve van Der Mije (Naturalis Biodiv Center, Leiden, Netherland), Burton Lim (Royal Ontario Museum Toronto, Canada), and François Catzeflis (UM-ISEM). We also thank Alison Devault (Arbor Biosciences) for her help on the capture kit design. We acknowledge the financial support of the "Agence Nationale de la Recherche" through the GAARanti project grant (ANR-17-CE31-0009). P.-H.F. acknowledges the SYNTHESYS Foundation for funding his work in the BMNH collections (GB-TAF-2735, GB-TAF-5026, and GB-TAF-5737). P.-H.F. also obtained financial support from Agence Nationale de la Recherche (CEBA, ref. ANR-10-LABX-25-01). YLRL had continuous support from Fundação de Amparo a Pesquisa e Inovação do Espírito Santo (FAPES) and Conselho Nacional de Desenvolvimento Científico e Tecnológico (CNPq), Brazil. The bio-informatics analyses of this article largely benefited from the Montpellier Bioinformatics Biodiversity cluster computing platform. This publication is the contribution n°XXX-XXX SUD of the Institut des Sciences de l'Évolution de Montpellier (UMR 5554 UM-CNRS-IRD-EPHE).

Authors declare no competing interests. TABLES

Species

GC3 content Relative evolution rate Exon coverage

Distance to reference