Metagenomics uncovers dietary adaptations for chitin digestion in the gut microbiota of convergent myrmecophagous mammals Sophie Teulleta,#, Marie-Ka Tilaka, Amandine Magdeleinea, Roxane Schaubb,c, Nora M. Weyerd, Wendy Panainod,e, Andrea Fullerd, William. J. Loughryf, Nico L. Avenantg, Benoit de Thoisyh,i, Guillaume Borrelj and Frédéric Delsuca,# aInstitut des Sciences de l’Evolution de Montpellier (ISEM), Univ Montpellier, CNRS, IRD, Montpellier, France bCIC AG/Inserm 1424, Centre Hospitalier de Cayenne Andrée Rosemon, Cayenne, French Guiana cTropical Biome and immunopathology, Université de Guyane, Labex CEBA, DFR Santé, Cayenne, French Guiana dBrain Function Research Group, School of Physiology, University of the Witwatersrand, Johannesburg, South Africa eCentre for African Ecology, School of Animals, Plant, and Environmental Sciences, University of the Witwatersrand, Johannesburg, South Africa fDepartment of Biology, Valdosta State University, Valdosta, GA, USA gNational Museum and Centre for Environmental Management, University of the Free State, Bloemfontein, South Africa hInstitut Pasteur de la Guyane, Cayenne, French Guiana, France iKwata NGO, Cayenne, French Guiana, France jInstitut Pasteur, Université Paris Cité, UMR CNRS 6047, Evolutionary Biology of the Microbial Cell, Paris, France #Corresponding authors: sophie.teullet@umontpellier.fr; frederic.delsuc@umontpellier.fr Abstract In mammals, myrmecophagy (ant and termite consumption) represents a striking example of dietary convergence. This trait evolved independently at least five times in placentals with myrmecophagous species comprising aardvarks, anteaters, some armadillos, pangolins, and aardwolves. The gut microbiome plays an important role in dietary adaptation, and previous analyses of 16S rRNA metabarcoding data have revealed convergence in the composition of the gut microbiota among some myrmecophagous species. However, the functions performed by these gut bacterial symbionts and their potential role in the digestion of prey chitinous exoskeletons remain open questions. Using long- and short-read sequencing of fecal samples, we generated 29 gut metagenomes from nine myrmecophagous and closely related insectivorous species sampled in French Guiana, South Africa, and the USA. From these, we reconstructed 314 high-quality bacterial genome bins of which 132 carried chitinase genes, highlighting their potential role in insect prey digestion. These chitinolytic bacteria belonged mainly to the family Lachnospiraceae, and some were likely convergently recruited in the different myrmecophagous species as they were detected in several host orders (i.e., Enterococcus faecalis, Blautia sp), suggesting that they could be directly involved in the adaptation to myrmecophagy. Others were found to be more host-specific, possibly reflecting phylogenetic constraints and environmental influences. Overall, our results highlight the potential role of the gut microbiome in chitin digestion in myrmecophagous mammals and provide the basis for future comparative studies performed at the mammalian scale to further unravel the mechanisms underlying the convergent adaptation to myrmecophagy. Main figures and corresponding datasets Figure_1_dataset.zip contains: FIGURE 1. Phylogenetic position of the 314 high-quality selected bins reconstructed from 29 gut metagenomes of the nine focal myrmecophagous species within a reference prokaryotic phylogeny. A: Phylogeny of the 314 selected bins (red branches) with 2496 prokaryote reference genomes. Circles respectively indicate (from inner to outer circles): the bacterial phyla and kingdom to which these genome bins were assigned based on the Genome Taxonomy Database release 7 (Parks et al, 2021). Clades, where a subtree was defined, are highlighted in blue for the Firmicutes (Fig. 1B), green for the Bacteroidetes, and pink for the Proteobacteria (Figs. S2 A and B, respectively). B: Subtree within Fimircutes showing myrmecophagous-specific clades (blue highlights; dark blue corresponds to the three clades mentioned in the results, light blue to the other clades). The outer circle indicates the bacterial family to which these genome bins were assigned based on the Genome Taxonomy Database. Bins’ names of the myrmecophagous-specific clades are indicated at leaves of the phylogenetic tree together with the genus to which they were assigned to. phylophlan_LR_SR_ToL_FINAL_concatenated.aln: Alignment of the concatenated markers assembled by PhyloPhlAn v3.0.58. phylophlan_LR_SR_ToL_FINAL.tre: Phylogenetic tree reconstructed by PhyloPhlAn v3.0.58 for the 314 high quality selected genome bins and the 2496 prokaryote reference genomes. Figure_2_dataset.zip contains: FIGURE 2. Phylogeny of the 394 GH18 sequences identified in 132 high-quality selected bins reconstructed from 29 gut metagenomes of the nine focal myrmecophagous species and relatives. Red branches indicate the 237 sequences having an active chitinolytic site (DXXDXDXE). Circles respectively indicate (from inner to outer circles): the bacterial family and phyla of the bin the sequence was retrieved from. Colored sequence names indicate the host species. Colored circles at certain nodes indicate enzymes to which sequences are similar when blasting them against the NCBI non-redundant protein database. Sequence names are indicated at leaves of the tree and begin with the genus to which the bin they were identified in was assigned to. GH18_sequences_from_selected_bins_alignment.fasta: Alignment of the 394 GH18 sequences identified in 132 high quality selected bins computed with MAFFT v7.450. GH18__sequences_from_selected_bins_tree.newick: Phylogenetic tree of the 394 GH18 sequences inferred with RAxML v8.2.11 within Geneious Prime 2022.0.2. Figure_3_dataset.zip contains: FIGURE 3. Detection of the 314 high-quality bacterial genomes (lines) in the 29 gut metagenomes (columns) of the nine focal species. Each square indicates the detection of a genome bin in a sample as estimated by anvi’o v7 (Eren et al, 2021). Names of bins are indicated on the left with red indicating chitinolytic bins (Table S2). The names begin with the genus to which the bin was assigned to. Asterisks (*) indicate bins detected in at least one soil sample (detection > 0.25) (Fig. S4, Table S2, and detection table available via Zenodo). Phylogenetic relationships of host species distinguished by different color strips are represented at the bottom of the graph. Columns on the right indicate (from left to right): the number of GH18 sequences identified in each bin (from 0 to 17), the bin’s taxonomic phylum, class, order, and family. The phylogeny of the 314 selected bins inferred with PhyloPhlAn v3.0.58 (Asnicar et al, 2020) is also represented on the right of the graph (see Fig. S1). Silhouettes were downloaded from phylopic.org. detection_bins_across_gut_metagenomes.txt: Detection table as tab-delimited file containing the detection values inferred by anvi'o v7 for the 314 high quality selected bins across the 29 gut metagenomes from the nine focal myrmecophagous species. Figure_4_dataset.zip contains: FIGURE 4. Distribution of chitinolytic selected bins (red links) among the nine focal myrmecophagous species and relatives. Phylogenies of the 314 high-quality selected bins (Fig. S1) and of the nine host species (downloaded from timetree.org) are represented respectively on the left and the right of the graph. Links illustrate, for each bin, in which host species the bin was detected (detection threshold > 0.25). Red links indicate bins in which at least one GH18 sequence with an active chitinolytic site (DXXDXDXE) was found (chitinolytic bins). The size of the circles at the tips of the host phylogeny is proportional to the number of samples (n = 1 for D. kap; n = 2 for D. nov, C. uni and M. tri; n = 3 for T. tet and O. afe; n = 4 for D. sp. nov FG; n = 6 for P. cri and S. tem). Bins’ names are indicated at the tip of the bins’ phylogeny and main bacterial phyla are indicated by colored vertical bars. This graph was done with the cophylo R package within the phytools suite (Revell, 2012). Silhouettes were downloaded from phylopic.org. presence_absence_MAGs_in_metagenomes.txt: Presence/absence matrix of the 314 selected genome bins across the 29 gut metagenomes. host_species_phylo_reduced_fig4.newick: Host phylogenetic timetree. Table_1_sample_infos.xls: Detailed sample information for the 33 fecal samples collected. N.B.: Diet was determined based on field observations (i.e., dissections) and the literature. Supplementary results Supplementary_results_Teullet_etal_2023.zip includes a comparison of genome statistics of the selected bins reconstructed from the long-read vs the short-read datasets, a phylogeny of the set of selected bins before dereplication (n = 407) and a comparison of the distribution of shared and specific genome bins carrying GH18 among host orders. Supplementary material Supplementary_material_Teullet_etal_2023.zip contains Supplementary figures (S1-S4) and tables (S1-S4). phylophlan_314_bins_phylogeny_FINAL_concatenated.aln and phylophlan_314_bins_phylogeny_FINAL.tre: Alignment of the concatenated markers and the final tree (respectively) reconstructed by PhyloPhlAn v3.0.58 for the 314 high-quality selected and dereplicated genome bins. phylophlan_407_selected_bins_nodRep_concatenated.aln and phylophlan_407_selected_bins_phylogeny_FINAL.tre: Alignment of the concatenated markers and the final tree (respectively) reconstructed by PhyloPhlAn v3.0.58 for the 407 high-quality selected genome bins before dereplication. abundance_bins_across_gut_metagenomes.txt: A tab-delimited file corresponding to the absolute abundance values inferred by anvi'o v7 for the 314 high-quality selected bins across the 29 gut metagenomes from the nine focal myrmecophagous species. detection_bins_across_soil_samples.txt: A tab-delimited file corresponding to the detection values inferred by anvi'o v7 for the 140 high-quality selected bins reconstructed from the aardvark, ground pangolin and southern aardwolf gut metagenomes across the eight soil samples collected on sample sites in South Africa. Assemblies Long-read_metagenomic_assemblies_polished.zip contains the 31 long-read metagenomes assembled with metaFlye strain v2.9 and polished with short reads using Pilon v1.4, which were used for binning. Long-read_metagenomic_assemblies_not_polished.zip contains the 33 long-read metagenomes assembled with metaFlye strain v2.9 before polishing. Short-read_metagenomic_assemblies.zip contains the 31 short-read metagenomes assembled with metaSPAdes and MEGAHIT. N.B: Two samples (DASY M1746 and DASY VLD168) were not sequenced using Illumina short reads. Only long reads were generated and assembled for these two samples and are made available here. As these assemblies could not be polished, these samples were not included in downstream analyses. Two samples (CAB M3141 and MYR M5293) were highly contaminated by host reads and not used in downstream analyses. As they were still assembled with the other samples, the corresponding metagenomes are made available here. Binning: genome bins and dereplication results High-quality_selected_bins_dereplicated.zip contains the 314 high quality selected bins (>90% completion, <5% redundancy) reconstructed from long- and short-read metagenomes with metaBAT2 and dereplicated with dRep at 98% ANI. metaBAT2_short-read_assemblies_bins.zip contains all bins reconstructed from the short-read assemblies with metaBAT2 (i.e., output of metaBAT2). metaBAT2_long-read_assemblies_bins.zip contains all bins reconstructed from the long-read polished assemblies with metaBAT2 (i.e., output of metaBAT2). Output_dRep_98ANI_407_bins_long-short-reads.zip contains the output of the dereplication analysis done on the set of 407 high-quality selected genome bins reconstructed from long- (n = 201) and short-read (n = 206; labeled "spad") metagenomes. It was performed with dRep using default parameters. After this step, the final dataset included 314 high-quality non-redundant genome bins. This folder includes: LR_SR_407_bins_dRep_98ANI_Primary_clustering_dendrogram.pdf: The primary clustering of selected genome bins using the Mash algorithm with an ANI threshold of 90%. LR_SR_407_bins_dRep_98ANI_Secondary_clustering_dendrograms.pdf: The secondary clustering of selected genome bins using the fastANI algorithm with an ANI threshold of 98%. LR_SR_407_bins_dRep_98ANI_Cluster_scoring.pdf: The clustering score attributed to each genome bin during dereplication. Asteriks (*) indicate genomes chosen to be the representative genomes of their cluster.