The nine-banded armadillo (Dasypus novemcinctus) is the most widespread xenarthran species across the Americas. Recent studies have suggested it is composed of 4 morphologically and genetically distinct lineages of uncertain taxonomic status. To address this issue, we used a museomic approach to sequence 80 complete mitogenomes and capture 997 nuclear loci for 71 Dasypus individuals sampled across the entire distribution. We carefully cleaned up potential genotyping errors and cross-contaminations that could blur species boundaries by mimicking gene flow. Our results unambiguously support 4 distinct lineages within the D. novemcinctus complex. We found cases of mito-nuclear phylogenetic discordance but only limited contemporary gene flow confined to the margins of the lineage distributions. All available evidence including the restricted gene flow, phylogenetic reconstructions based on both mitogenomes and nuclear loci, and phylogenetic delimitation methods consistently supported the 4 lineages within D. novemcinctus as 4 distinct species. Comparable genetic differentiation values to other recognized Dasypus species further reinforced their status as valid species. Considering congruent morphological results from previous studies, we provide an integrative taxonomic view to recognize 4 species within the D. novemcinctus complex: D. novemcinctus, D. fenestratus, D. mexicanus, and D. guianensis sp. nov., a new species endemic of the Guiana Shield that we describe here. The 2 available individuals of D. mazzai and D. sabanicola were consistently nested within D. novemcinctus lineage and their status remains to be assessed. The present work offers a case study illustrating the power of museomics to reveal cryptic species diversity within a widely distributed and emblematic species of mammals.

Here is a brief description of the files present in each folder:List of genes: List Immunity genes: Description and references for the immunity genes. list_Sma3s: List of genes identified by Sma3s program list_database: List of genes identified Immunome Knowledge Base, InnateDB and gene annotation. Polymorphism analysis : README: help to run seq_stat_coding seq_stat_coding: Executable C++ file and source file (.cpp) to estimate synonymous nucleotide genetic diversity (Ps) and non-synonymous nucleotide diversity (Pn) removeStopCodon: Executable C++ file and source file (.cpp) to remove the stop codon from the alignments Clean_Alignment : Executable C++ fileand source file (.cpp) to exclude the site with less than a define number of individuals TLR7_Pmajor_XP_015492521.1: Input file example Scripts and dataset : Simulations with SliM: SliM__immunity.slim: Script used to simulate sequences under balancing selection using codon format and exon-intron format Formating dataset: out_Seq_stat_*: Output of Seq_stat_coding according to our differents gene categories and selection pressure Formating_dataframe_script: Script to calculate Pn/Ps ratio and make a dataframe ready to plot. Plots and models : Dataframe_ready_to_plot.csv: Output from Formating_dataframe_script. Table containing : Species, PNPS, PN, PS, D_Taj, GC, S, category, Origin, family and selection regime. Plot_and_models_script: Script to plot results, and make differents models data_from_Leroy_etal_2021.csv : Informations and statistics about dataset from Leroy et al., 2021. Mitochondrial_phylogeny.treefile: Phylogeny based on mitochondrial genes of species from the dataset reconstructed by maximum likelihood method (IQTREE model GTR+Gamma and ultrafast bootstrap). Simulation analysis Tab_h_overdominance : Effect of parameters h (dominance coefficient) on Pn/Ps for sequences simulated under overdominace. Tab_m3_freq_dep : Effect of parameters S (selection coefficient) on Pn/Ps for sequences simulated under frequency dependent. Tab_Ne_overdominance : Effect of population size on Pn/Ps for sequences simulated under overdominance Tab_Ne_freq_dep : Effect of population size on Pn/Ps for sequences simulated under frequency dependent. Plot_simulated_results: Script to plot the effect of parameters on Pn/Ps from simulations. Supplementary table and figures : Figure S1 : Distribution of the percentage of contaminating contigs. The red line represents the 80% quantile. Figure S2: PCA of Cyanistes species Figure S3: PCA of Cyanomitra and Turdus species Figure S4: PCA of Ploceus species Figure S5: Fis ~Nucleotide diversity Figure S6: Correlation between Pn/Ps (a) and Ps (b) calculated on the control genes in this study's dataset and those calculated by Leroy et al. (2021). Figure S7 - Missing_data_ps_Phylloscopus.pdf : Relationship between the maximum number of missing individuals allowed and synonymous nucleotide diversity (Ps) Phylloscopus trochilus and Fringilla coelebs. Figure S8: Effect of sub-sampling size on PN/PS Figure S9: Pn/PS according to Ps for sub-sampling control, TLR and BD genes. Figure S10: Boxplot of Pn/Ps according to population size for simulated sequences under overdomiance via SLiM Figure S11: Boxplot of Pn/Ps according to population size for simulated sequences under frequency dependence via SLiM Figure S12: Boxplot of Pn/Ps according to a) initial selection coefficient of the mutation under frequency dependence b) dominance coefficient (h) for simulated sequences under overdominance via SLiM Table S1 - Model comparison using reduce number of families : Model selection of all genes categories using reduce number of families (we grouped Turdidae within Muscicapidae, Nectariniidae, and Estrildidae within Ploceidae and Fringillidae within Thraupidae). Table S2 - Lm & PGLS on dPn/Ps : Alternative models (Lm for linear models and PGLS for Phylogenetic Generalized Least Squares) using the difference between Pn/Ps of immunity genes and control genes (Pn/Ps) as dependent variable, and species origin as explanatory variable. Table S3 - Samples & sequencing information : Table with sampling and sequencing information regarding the samples newly-sequenced in this study, and those obtained from Leroy et al. 2021 Table S4 - Quality of sequences per individual : Table with sequencing quality information ( number of genes analysed, proportion of available positions, depth coverage) regarding the samples newly-sequenced in this study, and those obtained from Leroy et al. 2021 Table S5 to S14 : Model selection by AICc criterion and ANOVA test. Summary of the best models. Mitochondrial phylogeny: AllSp_ultrafastaboot.treefile: The species phylogeny was estimated using mitochondrial genes and a maximum likelihood method implemented in IQTREE (model GTR+Gamma and ultrafast bootstrap AllSp_mito.fst : Alignment of the mitochondrial genes in fasta format. Table and figure of the main text: Figure 1: Phylogeny based on mitochondrial genes of species from the dataset Figure 2 : Conceptual diagram showing the expected results Figure 3 : Boxplot of Pn/Ps according to species origin for different gene categories under purifying selection. Figure 4 : Boxplot of Pn/Ps according to species origin for different gene categories under purifying selection. Figure 5 : Effect of mutation type on Pn/Ps acording to Ne Figure 6 : Boxplot of Pn/Ps according to species origin (mainland in green and insular in orange) for different gene categories under balancing selection. Table 1: List of species and sampling localities, along with the type of data obtained and the number of individuals (N). Table 2 : Statistical model explaining Pn/Ps variation of Toll-Like Receptors, Beta-Defensins genes, and control genes. The p-values of ANOVA test between simpler models are not reported if a more complex model explains a larger proportion of the variance. Table 3 : Summary of the best statistical model selected using AICc explaining variation in Pn/Ps in control genes, Toll-Like receptors and Beta-Defensins genes under purifying selection with origin, gene category parameters.* indicates significances : * < 0.05; ** < 0.01; *** < 0.001. Table 4 : Statistical model explaining Pn/Ps variation of genes under balancing selection (i.e MHC class I and II), and simulated sequences under frequency dependence.

Shared ecological conditions encountered by species that colonize islands often lead to the evolution of convergent phenotypes, commonly referred to as “island syndrome”. Reduced immune functions have been previously proposed to be part of the island syndrome, as a consequence of the reduced diversity of pathogens on island ecosystems. According to this hypothesis, immune genes are expected to exhibit genomic signatures of relaxed selection pressure in island species. In this study, we used comparative genomic methods to study immune genes in island species (N = 20) and their mainland relatives (N = 14). We gathered public data as well as generated new data on innate (Toll-Like Receptors, Beta Defensins) and acquired immune genes (Major Histocompatibility Complex classes I and II), but also on hundreds of genes annotated as involved in various immune functions. As a control, we used a set of 97 genes not involved in immune functions, to account for the lower effective population sizes in island species. We used synonymous and non-synonymous variations to estimate the selection pressure acting on immune genes. For the genes evolving under balancing selection, we used simulation to estimate the impact of population size variation. We found a significant effect of drift on immune genes of island species leading to a reduction in genetic diversity and efficacy of selection. However, the intensity of relaxed selection was not significantly different from control genes, except for MHC class II genes. These genes exhibit a significantly higher level of non-synonymous loss of polymorphism than expected assuming only drift and an evolution under frequency dependent selection, possibly due to a reduction of extracellular parasite communities on islands. Overall, our results showed that demographic effects lead to a decrease in the immune functions of island species, but the relaxed selection caused by a reduced parasite pressure may only occur in some immune genes categories.