Genomic diversity and ecology of human-associated Akkermansia species in the gut microbiome revealed by extensive metagenomic assembly

Abstract

Background: Akkermansia muciniphila is a human gut microbe with a key role in the physiology of the intestinal mucus layer and reported associations with decreased body mass and increased gut barrier function and health. Despite its biomedical relevance, the genomic diversity of A. muciniphila remains understudied and that of closely related species, except for A. glycaniphila, unexplored.

Results: We present a large-scale population genomics analysis of the Akkermansia genus using 188 isolate genomes and 2226 genomes assembled from 18,600 metagenomes from humans and other animals. While we do not detect A. glycaniphila, the Akkermansia strains in the human gut can be grouped into five distinct candidate species, including A. muciniphila, that show remarkable whole-genome divergence despite surprisingly similar 16S rRNA gene sequences. These candidate species are likely human-specific, as they are detected in mice and non-human primates almost exclusively when kept in captivity. In humans, Akkermansia candidate species display ecological co-exclusion, diversified functional capabilities, and distinct patterns of associations with host body mass. Analysis of CRISPR-Cas loci reveals new variants and spacers targeting newly discovered putative bacteriophages. Remarkably, we observe an increased relative abundance of Akkermansia when cognate predicted bacteriophages are present, suggesting ecological interactions. A. muciniphila further exhibits subspecies-level genetic stratification with associated functional differences such as a putative exo/lipopolysaccharide operon.

Conclusions: We uncover a large phylogenetic and functional diversity of the Akkermansia genus in humans. This variability should be considered in the ongoing experimental and metagenomic efforts to characterize the health-associated properties of A. muciniphila and related bacteria.

Fig. 1

The Akkermansia genus comprises four additional candidate species phylogenetically rooted between the already characterized A. glycaniphila and A. muciniphila. A Whole-genome phylogeny of the 2420 metagenome-assembled genomes (MAGs) reconstructed here and the genomes from isolate sequencing available in NCBI taxonomically annotated as A. muciniphila or Akkermansia spp. The phylogenetic tree is rooted using Verrucomicrobium spinosum as an outgroup and was built using PhyloPhlAn 3 [46] with 400 universal markers (see the “Methods” section). SGB, species-level genome bin (see the “Methods” section). B Within- and between-clade whole-genome average estimated nucleotide identity (fastANI [47], top panels) and full-length 16S sequence distances (bottom panels) among Akkermansia SGBs provide evidence that these are candidate species

Fig. 2

Prevalence and insights into the ecological and functional characteristics of Akkermansia candidate species. A Akkermansia candidate species have variable prevalence across hosts and wild versus captive mice and non-human primates. We computed prevalences using species-specific marker genes (see the “Methods” section) applied on a total of 13,237 metagenomic samples. B, C Akkermansia candidate species are strongly mutually exclusive (analysis based on 4171 Akkermansia-positive human metagenomes). D A. muciniphila but not the other Akkermansia candidate species is associated with decreased host body mass index (BMI) according to a meta-analysis random effect model of partial correlations adjusted for age and sex (see the “Methods” section) comprising 3311 human metagenomic samples from 22 datasets (Additional file 1: Table S2). E Corrin ring biosynthesis operon genes are consistently present only in candidate species SGB9227 and SGB9228 (see the “Methods” section). F Growth analysis of the A. muciniphila and A. glycaniphila type strains shows propionate production by Pyt^T but not Muc^T in the absence of vitamin B12. This is indicative of endogenous production of vitamin B12 (acting as a cofactor for the methyl-malonyl CoA synthase reaction) by Pyt^T but not Muc^T. G Core gene genetic distances are correlated with corrin ring biosynthesis gene genetic distances. Pairwise distances were computed only for strains in which all genes were found together on the same contig

Fig. 3

The CRISPR-Cas system of Akkermansia candidate species and their viral targets. A CRISPR locus type composition of Akkermansia candidate species. All candidate species possess CRISPR locus type I-C, with the exception of A. muciniphila in which type II-C is present in more than 30% of the genomes. B Representative locus organization of CRISPR loci over Akkermansia candidate species. Some type I-C loci contain only one CRISPR array. Gene and CRISPR array lengths are scaled to correspond to the median length over all loci. C Phylogenetic tree of A. muciniphila subspecies colored by type II-C presence. D The total number of spacer sequences for the genomes in each Akkermansia candidate species. Type II-C loci were only found in A. muciniphila. Numbers above the boxplots correspond to the fraction of type I-C loci with two CRISPR arrays. E Logo plots of predicted PAM sequences in putative (phage) Viral Clusters (VCs, see the “Methods” section) upstream of sequences with perfect matches against CRISPR spacer sequences from type I-C loci. F Proportion of CRISPR spacers within candidate species genomes with a near-perfect match (at most 2 mismatched nucleotides) for four VCs. The number above the box plots corresponds to the fraction of genomes with at least one spacer hit against a given VC (see the “Methods” section). G Mapping of spacers from Akkermansia genomes against two representative VCs, visualized with a sliding window of 150 nt. See Additional file 2: Figure S8 for the remaining VCs. H Distribution of the relative abundances of the Akkermansia candidate species based on the presence or absence of each cognate VC in the metagenome (Additional file 1: Table S2, see the “Methods” section). P-values for differential abundance were determined via two-sided Wilcoxon rank-sum tests. P-values of <0.01 were considered significant. The numbers above the box plots correspond to the generalized fold change, with negative numbers indicating a higher bacterial abundance when a VC is detected [60]

Fig. 4

A. muciniphila is stratified in multiple subspecies with distinct host preferences. A Phylogenetic tree of A. muciniphila based on a core-gene alignment built using 169 clade-specific core genes (see the “Methods” section). The red arrow indicates the Muc^T type strain. B Within- and between-subspecies core-gene nucleotide identities confirm the subspecies diversification defined on the phylogeny. C Per-host frequency of A. muciniphila subspecies assembled from metagenomes. All 174 mouse A. muciniphila genomes were reconstructed from stool metagenomes of laboratory mice

Fig. 5

Functional diversification of A. muciniphila subspecies and cognate exopolysaccharide/LipidA synthesis operon. A Ordination analysis (Jaccard-distance-based PCoA using gene presence and absence information) reveals a diversification of gene repertoires of A. muciniphila subspecies. Genes found in less than 3% of strains were excluded. Subspecies designation is derived from the A. muciniphila phylogenetic tree in Fig. 4. B Operon archetypes putatively involved in exopolysaccharide/LipidA synthesis in A. muciniphila GP41 (operon archetype 1) and A. muciniphila Muc^T (operon archetype 2). C PCoA (same as in A) colored by operon archetype membership. Genomes in which neither operon could be found were labeled “Unassigned”