DOI: 10.5281/zenodo.4043409
title: High-quality carnivore genomes from roadkill samples enable species delimitation in aardwolf and bat-eared fox

description: <strong>High-quality carnivore genomes from roadkill samples enable species delimitation in aardwolf and bat-eared fox</strong> Rémi Allio<sup>1</sup>*, Marie-Ka Tilak<sup>1</sup>, Céline Scornavacca<sup>1</sup>, Nico L. Avenant<sup>2</sup>, Erwan Corre<sup>3</sup>, Benoit Nabholz1, and Frédéric Delsuc<sup>1</sup>* <strong>Affiliations</strong> <sup>1</sup>Institut des Sciences de l’Evolution de Montpellier (ISEM), CNRS, IRD, EPHE, Université de Montpellier, France remi.allio@umontpellier.fr marie-ka.tilak@umontpellier.fr celine.scornavacca@umontpellier.fr benoit.nabholz@umontpellier.fr frederic.delsuc@umontpellier.fr <sup>2</sup>National Museum and Centre for Environmental Management, University of the Free State, Bloemfontein, South Africa navenant@nasmus.co.za <sup>3</sup>CNRS, Sorbonne Université, FR2424, ABiMS, Station Biologique de Roscoff, 29680 Roscoff, France corre@sb-roscoff.fr *Correspondence: remi.allio@umontpellier.fr, frederic.delsuc@umontpellier.fr <strong><em>Running head</em></strong> Genomics from roadkill samples <strong>Abstract</strong> In a context of ongoing biodiversity erosion, obtaining genomic resources from wildlife is becoming essential for conservation. The thousands of yearly mammalian roadkill could potentially provide a useful source material for genomic surveys. To illustrate the potential of this underexploited resource, we used roadkill samples to sequence reference genomes and study the genomic diversity of the bat-eared fox (<em>Otocyon megalotis</em>) and the aardwolf (<em>Proteles cristata</em>) for which subspecies have been defined based on similar disjunct distributions in Eastern and Southern Africa. By developing an optimized DNA extraction protocol, we successfully obtained long reads using the Oxford Nanopore Technologies (ONT) MinION device. For the first time in mammals, we obtained two reference genomes with high contiguity and gene completeness by combining ONT long reads with Illumina short reads using hybrid assembly. Based on re-sequencing data from few other roakill samples, the comparison of the genetic differentiation between our two pairs of subspecies to that of pairs of well-defined species across Carnivora showed that the two subspecies of aardwolf might warrant species status (<em>P. cristata</em>and <em>P. septentrionalis</em>), whereas the two subspecies of bat-eared fox might not. Moreover, using these data, we conducted demographic analyses that revealed similar trajectories between Eastern and Southern populations of both species, suggesting that their population sizes have been shaped by similar environmental fluctuations. Finally, we obtained a well resolved genome-scale phylogeny for Carnivora with evidence for incomplete lineage sorting among the three main arctoid lineages. Overall, our cost-effective strategy opens the way for large-scale population genomic studies and phylogenomics of mammalian wildlife using roadkill. <strong><em>Figures & Tables</em></strong> <strong>Figure 1</strong>. Disjunct distributions of the aardwolf (<em>Proteles cristata</em>) and the bat-eared fox (<em>Otocyon megalotis</em>) in Eastern and Southern Africa. Within each species, two subspecies have been recognized based on their distributions and morphological differences (Clark, 2005; Koehler and Richardson, 1990). <strong>Figure 2</strong>. Representation of the mitochondrial genetic diversity within Carnivora with a) the mitogenomic phylogeny inferred from 142 complete Carnivora mitogenomes including those of the two populations of aardwolf (<em>Proteles cristata</em>) and bat-eared fox (<em>Otocyon megalotis</em>) and b) intraspecific (orange) and the interspecific (red) genetic diversities observed for the two mitochondrial markers COX1 and CYTB. <strong>Figure 3. </strong>Comparison of 503 mammalian genome assemblies from 12 taxonomic groups using bean plots of the a) number of scaffolds, and b) scaffold N50 values ranked by median values. Thick black lines show the medians, dashed black lines represent individual data points, and polygons represent the estimated density of the data. Note the log scale of the Y axes. The bat-eared fox (<em>Otocyon megalotis</em>) and aardwolf (<em>Proteles cristata</em>) assemblies produced in this study using SOAPdenovo and MaSuRCA are indicated by asterisks. Bean plots were computed using BoxPlotR (Spitzer et al., 2014). <strong>Figure 4. </strong>BUSCO completeness assessment of 67 Carnivora genome assemblies visualized as bar charts representing percentages of complete single-copy (light blue), complete duplicated (dark blue), fragmented (yellow), and missing (red) genes ordered by increasing percentage of total complete genes. The bat-eared fox (<em>Otocyon megalotis</em>) and aardwolf (<em>Proteles cristata</em>) assemblies produced in this study using MaSuRCA and SOAPdenovo are indicated by asterisks. <strong>Figure 5. </strong>Genetic differentiation indices obtained from the comparison of intraspecific (orange) and interspecific (red) polymorphisms in four pairs of well-defined Carnivora species and for the subspecies of aardwolf (<em>Proteles cristata</em>) and bat-eared fox (<em>Otocyon megalotis</em>) (grey). <strong>Figure 6. </strong>PSMC estimates of the change in effective population size over time for the Eastern (orange) and Southern (blue and purple) populations of a) bat-eared fox and ) aardwolf. mu = mutation rate of 10<sup>-8</sup>mutations per site per generation and g = generation time of 2 years. Vertical red lines indicate 20kyrs and 40kyrs. <strong>Figure 7. </strong>Phylogenomic tree reconstructed from the nucleotide supermatrix composed of 14,307 single-copy orthologous genes for 52 species of Carnivora plus one outgroup (<em>Manis javanica</em>). The family names in the legend are ordered as in the phylogeny. <strong>Table 1. </strong>Summary of sequencing and assembly statistics of the genomes generated in this study. <em><strong>Additional files</strong></em> <strong>Figure S1</strong>: Plot of the quality of Nanopore long reads base-called with either the <em>fast</em>or the <em>high accuracy</em>option of Guppy v3.1.5. The quality of the base-calling step has a large impact on the final quality of the assemblies by reducing the number of contigs and increasing the N50 value. <strong>Figure S2</strong>: Definition of the genetic differentiation index (GDI) based on the F-statistic (FST). The main difference between these two indexes is the use of heterozygous allele states for GDI rather than real polymorphism for the FST. Green = π<sub>within</sub>, Orange = π<sub>between</sub>, Blue = Population A, Red = Population A+B. <strong>Figure S3</strong>: Graphical representation of the results of contamination analyses performed with BlobTools for a) the aardwolf (<em>Proteles cristata</em>) and b) the bat-eared fox (<em>Otocyon megalotis</em>). <strong>Table S1</strong>: Pairwise patristic distances estimated for the 142 species based on the phylogenetic tree inferred with the 15 mitochondrial loci (2 rRNAs and 13 protein-coding genes). <strong>Table S2</strong>: Results of Bayesian dating for the two nodes leading to the <em>Proteles cristata</em>spp. and the<em>Otocyon megalotis</em>spp.. Divergence time estimates based on UGAM and LN models are reported with associated 95% credibility intervals for each MCMC chain. <strong>Table S3</strong>: Sample details and assembly statistics (Number of contigs/scaffolds and associated N50 values) for the 503 mammalian assemblies retrieved from NCBI (https://www.ncbi.nlm.nih.gov/assembly) on August 13th, 2019 with filters: “Exclude derived from surveillance project”, “Exclude anomalous”, “Exclude partial”, and using only the RefSeq assembly for <em>Homo sapiens</em>. <strong>Table S4</strong>: Genome completeness assessment of MaSuRCA and SOAPdenovo assemblies obtained for <em>Proteles cristata</em>and <em>Otocyon megalotis</em>together with the 63 carnivore assemblies available at NCBI and DNAZoo (https://www.dnazoo.org/assemblies) on August 13th, 2019 using Benchmarking Universal Single-Copy Orthologs (BUSCO) v3 with the Mammalia OrthoDB 9 BUSCO gene set. <strong>Table S5</strong>: Annotation summary and supermatrix composition statistics of the 53 species used to infer the genome-scale Carnivora phylogeny. <strong>Table S6</strong>:<strong> </strong>Sample details, Illumina sequencing, and assembly statistics of the 10 newly assembled carnivoran mitochondrial genomes. <strong>Table S7</strong>: Node calibrations used for the Bayesian dating inferences based on mitogenomic data. <strong>Table S8</strong>: Results of contamination analyses performed with BlobTools for the aardwolf (<em>Proteles cristata</em>). <strong>Table S9</strong>: Results of contamination analyses performed with BlobTools for the bat-eared fox (<em>Otocyon megalotis</em>). <strong>Table S10</strong>: Summary information for the Carnivora genomes available either on Genbank, DNAZoo (https://www.dnazoo.org) and the OrthoMaM database as of February 11th, 2020. The “OMM” column indicates if the genome was available on OMM (yes) or not (no). The “Annotation” column indicates whether the genome was already annotated (yes) or not (no). <strong><em>Zenodo supplementary files</em></strong> <strong>General_pipeline.sh</strong> contains the command lines used for the study 1- Mitogenomics Barcoding gap analyses: <strong>Barcoding gap analysis.zip</strong> contains the COX1 and CYTB matrices used for the barcoding gap analyses + the final trees and their associated IQ-TREE log files Mitochondrial phylogeny: <strong>Mitochondrial phylogeny.zip</strong> contains the mitochondrial supermatrix used to generate de phylogeny + the mitochondrial phylogeny and the associated IQ-TREE log file. Dating: <strong>Mitochondrial dating.zip</strong> contains the results of the dating analyses. 2- Genomics <strong>Assemblies.tar.gz</strong> contains the hybrid assemblies of <em>Proteles cristata cristata</em> and <em>Otocyon megalotis megalotis.</em> Busco analyses: <strong>Busco_Carnivora.xlsx</strong> contains genomes information like Accession numbers, Species names, Busco scores etc... Genetic differentiation: <strong>Genetic differentiation analyses.zip</strong> contains the coordinates of the regions used for the estimations for each genomes + a summary of the results for each species pair. PSMC analyses: <strong>PSMC.zip</strong> contains the results of the PSMC analyses. 3- Phylogenomics Phylogenomic analyses: <strong>Phylogenomic analysis.zip </strong>contains the nucleotides supermatrix composed of 14,307 single-copy orthologous genes for 54 species of carnivores plus one outgroup (<em>Manis javanica</em>) + the phylogenomic tree and the associated IQTREE log file. Gene trees: <strong>Gene trees.zip</strong> contains the gene trees.